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ABSTRACT 


The  buckling  load  of  a  blade  stiffened  laminated  composite  plate  having  midplane 
symmetry  is  maximized  for  a  given  total  weight.  The  thicknesses  of  the  layers  and  the 
width  and  height  of  the  stiffener  are  taken  as  the  design  variables.  Buckling  analysis  is 
carried  out  using  a  finite  element  method.  The  optimization  problem  is  solved  using 
commercially  available  optimization  packages. 

Due  to  the  highly  nonlinear  nature  of  the  optimality  equations,  several  local  opti¬ 
mum  solutions  are  found.  To  examine  the  relationship  between  the  number  of  local 
optimums  and  their  relative  magnitudes,  various  combinations  of  fiber  orientation  for 
the  laminate  layers  and  the  blade  stiffener  are  investigated. 
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The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made,  within 
the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic  er¬ 
rors,  they  cannot  be  considered  validated.  Any  application  of  these  programs  without 
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I.  INTRODUCTION 


Design  optimization  is  not  a  new  subject.  Man  has  always  strived  to  build  the  best, 
whether  it  be  ancient  man  making  a  hunting  bow  or  today's  engineer  designing  a  light¬ 
weight,  high-strength  wing  for  a  high-performance  aircraft.  The  goal  for  both  is  exactly 
the  same,  to  make  the  best  possible  design  subject  to  a  set  of  design  requirements.  The 
difference  between  the  two  is  not  in  the  goal  but  in  the  technique.  Ancient  man  used 
trial  and  error  to  improve  on  his  previous  designs  and,  over  a  period  of  time,  he  was  able 
to  produce  an  efficient  weapon.  Today,  high-speed  digital  computers  are  used  to  analyze 
and  improve  designs  by  solving  a  large  set  of  equations  that  have  been  formulated  to 
mathematically  model  the  design  problem. 

The  design  optimization  problem  of  current  interest  is  the  maximization  of  the 
buckling  load  for  a  blade  stiffened  composite  plate.  This  problem  can  be  viewed  from 
two  different  perspectives.  On  one  hand,  this  is  a  design  optimization  problem  which 
optimizes  the  shape  of  the  plate  and  stiffener  cross  section  to  give  the  maximum  load. 
On  the  other  hand,  this  is  a  design  optimization  problem  which  optimizes  the  thickness 
and/or  the  ply  angle  of  the  composite  laminae  to  give  the  maximum  load.  Previous 
works  on  the  optimization  of  the  plate  and  stiffener  cross  sectional  shape  have  focused 
on  the  use  of  solid  elastic  materials  (Refs.  1-4).  In  the  area  of  design  optimization  of 
composite  plates,  previous  works  have  focused  on  plates  without  stiffeners.  Optimiza¬ 
tion  of  the  composite  plate  with  respect  to  the  lamina  fiber  orientations  has  been  studied 
in  (Refs.  5-17).  Of  greater  interest  to  this  study  are  the  previous  wrorks  on  the  design 
optimization  of  composite  plates  where  the  laminae  thickness  and  not  the  ply  orien¬ 
tations  are  taken  as  the  design  variables  (Refs.  18-21).  The  current  study  combines  the 
design  optimization  problems  of  cross  section  optimization  and  composite  laminae 
thickness  optimization  into  a  single  design  optimization  problem.  This  is  done  by  opti¬ 
mizing  a  blade  stiffened  composite  plate  for  the  maximum  buckling  load  by  using  the 
laminae  thicknesses  and  the  stiffener  height  and  width  as  the  design  variables. 

Composites  are  used  as  the  structural  material  in  this  study  because  of  their  high 
strength-to-weight  and  stiffness-to-weight  ratios-  These  properties  are  very  important 
in  many  of  today's  engineered  structures,  especially  in  the  aerospace  industry  where  high 
strength  is  required  but  a  severe  penalty  is  incurred  for  added  weight.  Because  of  this, 
composites  make  an  excellent  choice  for  a  structural  material.  However,  the  use  of 
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composites  greatly  complicates  the  design  analysis  problem  arising  from  the  use  of  lam¬ 
inated  materials  with  non-homogeneous  properties.  Now  the  analysis  process  must 
consider  the  orthotropic  properties  and  the  ply  orientation  of  each  individual  lamina,  the 
specific  stacking  order  of  the  laminae,  and  the  thickness  of  each  individual  lamina.  One 
serious  problem  that  arises  is  the  coupling  effect  between  extension  and  bending.  Be¬ 
cause  a  symmetric  composite  plate  is  chosen,  the  coupling  between  extension  and 
bending  is  eliminated.  This  lack  of  coupling  allows  the  use  of  linear  buckling  analysis 
up  to  the  buckling  load. 

The  formulation  of  a  suitable  design  analysis  technique  is  the  first  step  in  the  design 
optimization.  Since  the  optimum  design  can  never  be  better  than  the  design  analysis 
that  it  is  based  on,  the  formulation  is  described  in  some  detail.  The  buckling  analysis 
of  the  stiffened  plate  is  performed  in  two  steps.  First,  a  finite  element  method  is  used 
to  determine  the  overall  buckling  load  of  the  stiffened  plate.  This  is  done  by  treating  the 
stiffener  as  a  beam  in  forming  the  global  stiffness  matrix  for  the  plate.  Using  FEM 
formulation,  the  first  three  buckling  loads  are  solved  numerically.  The  second  part  of 
the  buckling  analysis  finds  the  local  buckling  load  of  the  stiffener.  By  treating  the 
stiffener  as  a  plate  with  three  sides  simply  supported  and  the  fourth  side  as  free,  the  local 
buckling  problem  can  be  solved  analytically  using  a  Levy's  solution. 

Once  the  design  analysis  technique  has  been  formulated,  the  design  variables,  taken 
to  be  the  laminae  thickness  and  the  stiffener  height  and  width,  are  optimized  using  an 
optimization  program.  This  study  uses  three  different  optimization  programs  which  are 
all  commercially  available.  They  are:  the  1MSL  Subroutine  DNCONG;  the  Design 
Optimization  Tools  (DOT)  program  using  Modified  Feasible  Direction;  and  the  Design 
Optimization  Tools  (DOT)  program  using  Sequential  Linear  Programming. 

The  specific  problem  in  .stigated  in  this  study  is  a  symmetric  blade  stiffened  com¬ 
posite  plate  that  has  two  laminae  and  a  stiffener  on  either  side  of  the  centerline.  The 
fiber  orientation  of  the  laminae  and  the  stiffener  are  allowed  to  be  either  0°  or  90°.  Four 
different  stacking  orders  of  these  fiber  orientations  are  investigated.  For  each  case,  de¬ 
sign  optima  are  investigated.  These  optima  are  then  compared  to  determine  if  there  is 
a  best  stacking  configuration.  The  study  is  completed  for  two  different  total  volumes 
of  material.  Finally,  the  relative  efficiencies  of  the  three  optimization  routines  are  in¬ 
vestigated  to  determine  qualitatively  which  is  the  best  suited  for  this  specific  problem. 
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II.  DESIGN  APPROACH 


A.  GEOMETRY 

The  geometry  used  in  this  study  is  a  blade  stiffened  composite  plate  which  is  sym¬ 
metric  to  the  plate  centerline  as  shown  in  Figure  1  on  page  4.  The  plate  is  simply  sup¬ 
ported  on  all  four  sides  and  subject  to  a  compressive  inplane  load  in  the  stiffener 
direction.  The  blade  stiffener  is  centered  on  the  plate  and  spans  the  plate  length,  La. 

Each  lamina  has  a  thickness  of  T,  and  a  ply  orientation  of  angle  6,  .  The  stiffener 
is  defined  by  a  width  of  B  and  a  height  of  H.  The  ply  orientation  of  the  stiffener  can 
be  either  CF  or  90°  .  The  independent  engineering  properties  of  each  lamina  and  the 
stiffener  are  given  by  £,, ,  E^  ,  v,v  Gn  ,  and  E„4,  v,v  Gl2s  ,  respectively. 

B.  STRESS  ANALYSIS 

The  analysis  of  a  composite  plate  structure  is  based  on  classical  lamination  theory. 
From  this  theory,  the  stiffness  of  a  composite  structure  can  be  determined  from  the 
material  properties  of  the  individual  laminae  [Ref.  22:  p.  147).  The  plate  in  this  study  is 
composed  of  individual  layers  which  are  made  of  orthotropic  material.  Also,  the  plate 
is  thin  with  respect  to  its  span  length  so  it  is  assumed  to  be  under  plane  stress  conditions. 
From  this,  classical  lamination  theory'  can  be  applied  to  determine  the  composite's 
stiffness. 

First,  the  stress  strain  relation  for  a  orthotropic  material  under  plane  stress  must  be 
defined  for  the  individual  laminae.  This  relationship  is  given  in  Eq.  (1)  below. 


Q\\  Q\l  Q 16 
Q\2  0.22  C?26 
Q 16  C?26  Q(>(> 


(1) 


where  the  reduced  stiffness  matrix  elements,  Q,y,  are  defined  in  terms  of  the  elastic 
modulus  in  the  1-1  principal  direction,  £j,,  the  elastic  modulus  in  the  2-2  principal  di¬ 
rection,  Eu,  the  shear  modulus,  Gu ,  and  poisson  s  ratio  relating  transverse  strain  in  the 
1  direction  when  stressed  in  the  2  direction,  v12  ,  and  poisson  s  ratio  relating  transverse 
strain  in  the  2  direction  when  stressed  in  the  1  direction,  v21  . 
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Figure  1.  Problem  Geometry:  top)  general  view  bottom)  cross  sectional 

view 
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For  each  lamina,  the  principal  axes  lie  along  the  direction  of  the  fiber  orientation. 
Since  the  fiber  orientation  of  each  lamina  is  allowed  to  vary  by  an  angle  0,  Eq.  (1)  must 
be  transformed  to  determine  the  stress  strain  relationship  in  the  lamina  geometric  refer¬ 
ence  frame.  This  relationship  is  given  in  Eq.  (2). 


(2) 


where  the  elements  of  the  reduced  transformed  stiffness  matrix,  Q „  ,  are  defined  in  terms 
of  0V  and  9  as 

0n  =  0U  cos4®  +  2(012  +  2066)  sin2#  cos2#  +  022  sin4® 

012  =  (0u  +  022  ~  4066)  sin2®  cos2®  +  012(sin4®  +  cos4®) 

022  =  0ii  sin4®  +  2(0,2  +  2066)  sin2®  cos2®  +  022  cos4® 

016  =  (011  -  012  -  2066)  sin  9  cos3®  +  (0,2  -  022  +  2066)  sin3®  cos  ® 

026  “  (011  -  012  -  2066)  sin3®  cos  0  +  (012  -  022  +  2066)  sin  0  cos3® 

066  =  (011  +  022  -  20,2  -  2066}  sin3®  cos  0  +  066(  sin4®  +  cos4®) 


Now  Eq.  (2)  can  be  thought  of  as  the  stress  strain  relationship  for  the  A*  layer  of  a 
multilayered  laminate.  Thus  it  is  now  written  as  Eq.  (3). 
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M*  =  [QWk 


(3) 


The  strain  at  any  point  through  the  depth  of  the  laminate  can  be  expressed  in  terms  of 
the  midplane  strain,  e°,  the  distance  z  from  the  centerline,  and  the  midplane  curvature, 
k  .  See  Figure  2  on  page  7  for  the  definition  of  z  . 

The  strain  is  now  expressed  as  Eq.  (4). 


The  resultant  forces  and  moments  acting  on  a  laminate  are  obtained  by  integration 
of  the  stresses  in  each  layer  through  the  laminate  thickness. 
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Figure  2.  Half  of  Symmetric  Laminate  Plate 


Substituting  Eq.  (4)  into  Eqs.  (5a)  and  (5b)  gives, 
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_^16  ^26  ^66 _ 
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A 6  ^26  &66 

where 


n 

Aij  =  2X®*(Z*  ~  zk- 1) 


*=i 


n 

b,j  -  JjSi),  (*l  -  £,) 


k=\ 


(6) 


* 

A,  =  yX®*(z*  _  **-«) 

A=1 

The  A,,  terms  are  the  extensional  stiffness,  the  Bt/  terms  are  the  coupling  stiffness,  and 
the  DtJ  terms  are  the  bending  stiffness.  Since  the  composite  plate  modeled  in  the  study 
is  symmetric  about  the  middle  surface,  the  terms  B{J  all  drop  out  [Ref.  22:  p.  164], 
Therefore,  in  this  problem  there  is  no  coupling  between  extension  and  bending. 

C.  BUCKLING  ANALYSIS 

Buckling  of  the  stiffened  plate  can  occur  in  two  urays.  These  are:  the  plate  and 
stiffener  can  buckle  as  a  unit  into  one  of  the  global  buckling  modes;  and  the  blade 
stiffener  can  experience  local  buckling.  Because  of  this,  the  buckling  analysis  is  per¬ 
formed  in  two  steps.  For  the  first  case,  the  blade  stiffener  is  treated  as  a  beam  and  in¬ 
cluded  in  the  formulation  of  the  global  stiffness  matrix  in  the  Finite  Element 
Formulation.  For  the  second  case,  the  blade  stiffener  is  treated  as  a  plate  subject  to  a 
inplane  compressive  load  that  has  three  sides  simply  supported  and  the  fourth  side  free. 
[Ref.  23:  pp.  2-3J 

1.  Overall  Plate  Buckling  Analysis 

From  [Ref.  23:  p.  4],  the  overall  buckling  equation  in  matrix  form  is 

UQ{U)  +  [*]*{</}  ~  nJLKaMU)  -  pbl ATC],{U}  =  0 

where  [/Q  and  [K]>  are  the  plate  and  beam  stiffness  matrices,  respectively,  [Kc]  and 
[Kc]»  are  the  plate  and  the  beam  geometric  stiffness  matrices,  respectively,  and  {L’}  is  the 


8 


r'hl 


buckling  shape.  The  amount  of  load  carried  by  the  beam  and  plate  are  pb  and  n„  re¬ 
spectively.  The  load  share  carried  by  the  beam  and  by  the  plate  is  proportional  to  the 
relative  stiffness  of  each.  These  relative  loads  are  given  by 


epsp 


Vp  +  ebsb 


P 


Pb  = 


e^b 

epsp  +  ebsb 


P 


where  sp  and  sb  are  the  nondimensional  cross-sectional  area  of  the  plate  and  beam,  re¬ 
spectively,  and  ep  and  eb  are  the  nondimensional  elastic  moduli  of  the  laminated  plate 
and  the  beam,  respectively.  The  nondimensional  elastic  modulus  of  the  plate  must  be 
determined  from  the  properties  of  the  individual  laminae.  This  value  is  calculated  from 
the  extensional  stifTness  matrix  and  is  given  by  (Ref.  23:  p.  4]  as 

I  Ml 

(a22a66  -  al(>)lT 

Now  the  overall  buckling  equation  can  be  written  in  the  form, 

UQjiU)  -  plKaUU)  -  0 

where  [ K]T  is  the  total  stiffness  of  the  plate  and  stiffener  combination  and  [Afc]r  is  the 
total  geometric  stiffness  of  the  plate  and  stiffener  combination.  This  equation  can  easily 
be  recognized  as  an  eigenvalue  problem  with  p  representing  the  eigenvalues  and  {U}  re¬ 
presenting  the  eigenvectors.  This  problem  is  solved  using  DNLASO,  one  of  the  sub¬ 
routines  from  the  package  LAS02  (Ref.  24],  which  computes  a  few  eigenvalues  and  the 
associated  eigenvectors  of  a  large  (sparse)  symmetric  matrix  using  the  Lanczos  algorithm 
[Ref.  25].  The  programming  code  used  to  calculate  the  first  three  buckling  loads  of  the 
stiffened  plate  is  Subroutine  EIGEN  contained  in  Appendix  A. 
a.  FEM  Formulation 

In  order  to  determine  the  stiffness  matrices  for  the  buckling  equation,  a  fi¬ 
nite  element  formulation  is  required.  For  this  study,  a  16  degree  of  freedom  plate  ele¬ 
ment  is  chosen  for  the  plate  and  a  4  degree  of  freedom  beam  element  is  chosen  for  the 
stiffener.  The  plate  is  divided  into  four  elements  (2  X  2)  and  the  stiffener  is  divided  into 
two  elements. 
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From  (Ref.  26:  p.  i  I6j  the  beam  element  stiffness  matrix,  [/]4,  is  defined  as 


12 

-64 

-12 

-64" 

f»4 

-64 

442 

64 

2/2 

4J 

-12 

64 

12 

64 

_~64 

242 

64 

442_ 

where  the  element  length,  /„  is  one  half  the  nondimensionalized  span  length,  /„,  and  e„  is 
the  nondimensionalized  elastic  modulus  of  the  stiffener.  The  nondimensionalized  mo¬ 
ment  of  inertia  of  the  stiffener,  4,  is  defined  as 

4  =  y  b[{h  +  tT)3  -  i3T] 

where  tT  is  one  hall  ne  total  plate  thickness  for  the  symmetric  plate.  From  (Ref.  26:  p. 
388]  the  beam  element  geometric  stiffness  matrix,  [Ac]6,  is  defined  as 


~  36 

-34 

-36 

-34" 

Pb 

-34 

442 

34 
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34 
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34 
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For  the  plate  element,  the  geometric  stiffness  matrix  elements.  ,  are 
given  by  [Ref.  27:  p.  432]  as 


where  /  are  the  16  shape  functions  for  the  element.  The  shape  functions,/,  are  given 
below. 

/,  =  (1  -  3$2  +  2£3)(1  -  M2  +  2 r,3) 

h  =  (1  -  3S2  +  2£3)(>?  -  2V2  +  r,3) 

h  =  -  («  ”  2^2  +  $3)(l  -  3^2  +  2*3) 

h  -  (4  "  2^2  +  £3)(>?  -  2>?2  +  ^3) 
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A  -  -  2^3)(l  -  W  +  2^3) 


/«  =  (3£2  ~  2<^3) (?7  -  2t\2  +  >73) 

/7  =  (;2  -  ,‘3)(l  -  3,,2  +  2^) 

A  =  -  (;2  -  «%  -  2„2  +  O 

/9  =  (3^2  -  2{3)(3i,2  -  2„3) 
f\o  -  ~  (342  -  2S'3)(„2  -  „ 3) 
fn  =  (£2  -  «3)(3»f2  ~  2^3) 

fn  =  ~  W  -  >f3) 

fn  =  { I  -  3s'2  +  2c3)(3^2  -  2,3) 

/,4  =  ~(1  -  3?2  +  2^3)(^2  -  n') 

f\ 5  =  -(f  -  2;2  +  ^3)(3^2  -  2>/3) 

/■a  =  -  (4  -  2?2  +  ^  -  S) 

The  integral  of  Eq.  (7)  is  solved  numerically  by  using  gauss  quadrature. 

For  the  plate  element,  the  stiffness  matrix  elements,  k,n  are  given  by  [Ref. 
27:  pp.  415-418]  as 

kjj  —  f  I*  {^n^i  +  DnC2  +  £>16C3  +  D22C4  +  D26C5  +  D66C6}d^dn 
Jo  Jo 

where 
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Again,  here  the  integration  is  done  numerically  by  gauss  quadrature.  All  of  the  ele¬ 
mental  stiffness  formulation  is  done  in  Subroutine  ELESTF  in  Appendix  A. 

Now  that  all  the  elemental  stiffness  matrices  (both  stiffness  and  geometric) 
have  been  determined,  the  individual  elemental  matrices  can  be  assembled  into  the  global 
stiffness  matrices.  The  geometric  elemental  stiffness  matrices  [£c]  and  [/cc]4  are  com¬ 
bined  to  form  the  total  global  geometric  stiffness  matrix  [Kc]r .  The  elemental  stiffness 
matrices  [/c]  and  [£]*  are  combined  to  form  the  total  global  stiffness  matrix  [AT],.  This 
global  stiffness  matrix  assembling  is  done  in  Subroutine  ASSEMB  in  Appendix  A. 

2.  Stiffener  Local  Buckling  Analysis 

To  calculate  the  local  buckling  load  of  the  blade  stiffener,  it  is  treated  as  a  sep¬ 
arate  plate  buckling  problem.  Now  the  stiffener  alone  is  treated  as  a  plate  with  three 
sides  (x  =  0,  x  =  /„  z  =  0)  simply  supported  and  the  fourth  side  (z  =  h )  free.  It  is 
subject  to  axial  compressive  load  py  Because  of  the  geometry  and  the  fact  that  only  the 
first  buckling  lc  id  is  of  interest,  this  problem  can  be  solved  analytically. 
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The  general  partial  differential  equation  for  the  plate  buckling  problem  is  given 
from  (Ref.  22:  p.  260]  as 


Du  4t  +  +  2D, 

cx 


->4 

C  W 

6b/  3  2,  2 
cjr  cz 


+  D 


■3*4 

C  W 


22  ,  4 

CZ 


~  P  ' 


<3jc2 


=  0 


where  now  the  D  .  refer  to  the  properties  of  the  blade  stiffener  and  p  refers  to  the  local 
buckling  load  p„.  By  applying  the  boundary  conditions  and  assuming  a  Levy’s  solution 
of  the  form 


w 


=  V 


t  mux 

/  ,  zm  sin  — — 
*/? 


m=\  ,3,5.. 


the  buckling  problem  can  be  solved  directly  [Ref.  28:  p.  208]. 
earned  out  in  Subroutine  LEVYP  contained  in  Appendix  A. 


The  calculation  of  p ,  is 
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III.  OPTIMIZATION  PROBLEM 


A.  PROBLEM  STATEMENT 

The  optimization  problem  is  to  maximize  the  buckling  load  of  the  blade  stiffened 
piate  for  a  given  total  matenal  volume  (which  is  related  to  the  total  weight).  The  design 
variables  are  set  as  the  nondimensionalized  thickness  of  the  individual  lamina,  t„  and  the 
nondimensionalized  width  and  height  of  the  stiffener.  Here  the  nondimensionalized 
width,  b,  and  the  nondimensionalized  height,  h,  of  the  stiffener  are  denoted  as  and 
A..2  *  respectively.  The  nondimensional  design  variables,  t„  are  subject  to  side  constraints 
of 

'«mm  <  t,  <  U  mlx  for  i  =  1,2 . .(«  +  2) 

where  i,m„  and  i,min  are  upper  and  lower  bounds  on  the  design  variable  dimensions,  re¬ 
spectively,  and  n  is  half  the  number  of  layers  for  the  symmetric  plate. 

The  optimization  problem  for  maximizing  the  buckling  load  is  written  as 

max 

P 

b 

subject  to  2 £cr,  +  2^,,^  -0  =  0 

1.000/7,  >  p 
0.999/7,  >  p 

O. 998/7,  >  p 

P, >P 

and  t,  mm  <  r,  <r,  for  /=  1,2 (n  +  21 

where  p„  p2  and  p,  are  the  first  three  overall  buckling  loads  of  the  stiffened  plate,  pb  is  the 
local  buckling  load  of  the  stiffener,  t,  are  the  design  variables,  0  is  the 
nondimensionalized  total  plate  volume,  and  c  is  the  relative  w'idth  of  the  plate  LJL„.  P 
is  a  parameter  introduced  to  raise  all  the  buckling  loads  during  the  optimization  process, 
which  finally  converges  to  the  buckling  load  (lowest  eigenvalue)  of  the  optimized  design. 

The  first  constraint  represents  the  total  volume  constraint.  Simply  put,  the  total 
volume  of  the  design  cannot  exceed  the  maximum  available  resource.  The  next  four 
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constraints  are  minimum  limits  on  the  buckling  loads.  Here  the  buckling  load  of  the 
stiffened  plate  must  always  be  at  least  as  large  as  the  smallest  value  of  p„  p2  ,  p}  or  p 
The  coefficients  of  1.000,  0.999,  and  0.998  in  the  second,  third,  and  fourth  constraints, 
respectively,  are  necessary  to  allow  the  calculation  of  the  eigenvalues,  p{  ,  p2 ,  and  p2, 
when  the  buckling  mode  is  bimodal  or  trimoda!  [Ref.  29:  p.  355).  Finally,  the  side  con¬ 
straints  place  minimum  and  maximum  limits  on  the  size  of  each  design  variable. 

B.  OPTIMIZATION  METHODS 

The  optimization  problem  above  can  be  solved  by  a  number  of  different  methods. 
For  this  study  three  different  commercially  available  programs  are  used  to  solve  the  op¬ 
timization  problem.  Each  of  these  methods  uses  a  different  optimization  technique 
thereby  providing  three  independent  solutions  of  the  optimization  problem.  The  three 
methods  used  are:  the  IMSL  library  Subroutine  DNCONG;  Design  Optimization  Tools 
(DOT)  using  Modified  Feasible  Direction;  and  Design  Optimization  Tools  (DOT)  using 
Sequential  Linear  Programming.  Each  of  these  methods  are  variations  on  the  general 
optimization  routine.  The  general  steps  in  the  formulation  of  an  optimization  routine 
are  listed  below. 

1.  Stan  with  initial  guess  X°  for  the  0rt  iteration,  q  =  0. 

2.  Update  the  iteration  number,  q  =  q  +  1. 

3.  Evaluate  the  objective  function,  F(X),  and  the  constraint  functions,  g,(Y),  at  the 
current  value  of  X. 

4.  Identify  the  critical  constraints,  J. 

5.  Calculate  the  objective  function  gradient,  VFfiY),  and  the  critical  constraint  func¬ 
tion  gradients,  Vgy(X),  for  J  e  J. 

6.  Determine  the  search  direction,  S’. 

7.  Perform  a  one-dimensional  search  to  find  the  step  length,  a. 

8.  Update  the  solution  by  adding  the  previous  iterate  solution  to  the  product  of  the 
step  length  and  search  direction,  X*  —  X'~'  +  a'S’  . 

9.  Check  for  convergence.  YES:  done;  No:  go  to  step  2. 

For  each  of  the  optimization  routines  used  in  this  study,  a  description  of  the  variations 
from  the  standard  optimization  routine  are  given  below. 

1.  IMSL  Subroutine  DNCONG 

DNCONG  is  based  on  Subroutine  NLPQL  which  is  an  optimization  program 
developed  by  Schittkowski.  This  programming  algorithm  uses  a  successive  quadratic 
programming  method  to  solve  the  general  nonlinear  programming  problem.  In  this 
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method,  the  search  direction  subproblem  is  formulated  and  solved  by  using  a  quadratic 
approximation  of  the  Lagrangian  function  and  by  linearizing  the  constraints.  The 
Lagrangian  function,  L(X,  A),  is  given  by 


L(X,  /)  =  FIX'  -  Y^X) 
j=  i 

where  /  is  the  vector  of  Lagrange  multipliers  and  m‘  is  the  number  of  function  con¬ 
straints  plus  side  constraints,  m '  =  m  +  2n  .  The  search  direction  subproblem  is  for¬ 
mulated  as 

mn  n  ~-SrBqS  +  VftXfS 
SeUn  2  > 

subject  to  Vg/^)rS  +  gj(Xq)  =  0,  j  =  l,...,me 

VgfjFfS  +  gjiX*)  >  0,  1 

and  Xmn-X'  <  S  <  Xmix  -  X*  for  A  =1,2 . n 


where  B"  is  a  positive  definite  approximation  of  the  Hessian  matrix  of  the  Lagrangian 
function  and  X » is  the  current  iterate.  The  solution  of  the  subproblem,  S",  is  used  as  the 
search  direction  in  the  line  search  to  find  the  new  point  X ♦*l.  [Ref.  30) 

2.  DOT  using  MFD 

Design  Optimization  Tools  (DOT)  by  VMA  Engineering  is  a  programming  code 
designed  to  solve  a  variety  of  nonlinear  constrained  or  unconstrained  optimization 
problems.  The  version  of  DOT  that  uses  a  Modified  Feasible  Direction  algorithm  has 
a  search  direction,  S’,  that  is  a  modified  form  of  the  steepest  descent  direction.  The  ac¬ 
tual  modification  to  the  search  direction  is  done  by  using  the  Fletcher- Reeves  conjugate 
direction  method.  With  this  method,  the  search  direction  is  defined  as 

Sq  =  -\7F(Xq~')  +  C  Sq~' 


where 


c  = 


l7f(-OI 

Iw^r7-2)! 


2 


2 
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This  is  a  first  order  search  direction  which  is  extremely  simple  to  calculate.  The  modifi¬ 
cation  of  the  search  direction  gives  dramatically  improved  results  over  the  method  of 
steepest  descent.  [Ref.  31:  p.  E-16) 

3.  DOT  using  SLP 

Design  Optimization  Tools  (DOT)  using  Sequential  Linear  Programming  is  a 
method  that  approximates  a  nonlinear  problem  as  a  linear  problem  and  then  optimizes. 
First,  a  first  order  Taylor  Series  approximation  of  the  objective  and  constraint  functions 
are  calculated.  Then,  this  approximation  is  optimized  instead  of  the  original  nonlinear 
functions.  Since  the  problem  is  now  linear,  the  value  of  the  objective  and  constraint 
functions  are  easily  and  inexpensively  calculated.  Also,  now  the  gradients  of  the  objec¬ 
tive  and  constraint  functions  are  available  directly  from  the  Taylor  Series  expression. 
The  linear  approximation  problem  is  optimized  using  the  method  of  Modified  Feasible 
Direction.  [Ref.  31:  p.  E-33| 


IV.  RESULTS  AND  DISCUSSION 


The  laminate  plate  chosen  in  this  study  has  two  layers  with  equal  span  length,  /... 
and  span  width,  Lb.  The  material  selected  for  the  laminae  and  the  stiffener  is  a 
graphite  epoxy  composite.  The  material  properties  of  this  composite  are: 

=  31.0  x  10°  psi  ,  Eu  =  3.4  x  106  psi 

(7,2  =  0.75  x  106  psi  ,  v , 2  =  0.28 

Design  optimization  results  are  obtained  for  two  different  total  volume  conditions  (.0 
=  0.04  and  0.015)  and  four  different  fiber  orientation  stacking  sequences: 

( 0°beaml9Q°IO°)sym 
(0°beamIO°l90°)sym 
(90%,am/90°/0°)iym 
(90%„m/0 °l  90%m 

The  minimum  and  maximum  thickness  of  each  laminae  are  limited  to  0.00 IT,  and 
0.2  La,  respectively.  The  same  upper  and  lower  limits  are  used  for  the  stiffener  width  and 
height. 

A.  OPTIMUM  SOLUTIONS 

For  each  configuration,  the  known  optimums  and  the  values  of  the  design  variables 
at  these  optimums  are  given.  Also  listed  are  the  active  buckling  modes  and  the  number 
of  the  corresponding  figure  which  presents  the  same  information  graphically.  All  of  the 
results  are  presented  in  non-dimensional  form. 
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1.  Case  1:  8  =  0.04  and  Ply  Orientation  {0°^ml  90°/  0°),,,. 

For  this  configuration,  two  different  optimum  solutions  are  identified.  They  are 
listed  in  Table  1  below.  The  same  results  are  presented  graphically  in  Figures  3-6. 


Table  1.  RESULTS:  0  =  0.04,  PLY  ANGLE  (0 °UAMI  90°/  0°)sw 


<2 

b 

h 

P 

Buckling 

Modes 

0.00738 

0.01020 

0.03801 

0.06369 

0.0002646 

1.2. 3.4 

0.01521 

0.00100 

0.05794 

0.06526 

0.0002822 

1,2,4 

Figure  3.  Geometry  for  Case  1  with  p  =  0.0002646:  For  Case  1,  total  volume  0 
=  0.04,  with  ply  angle  orientation  (0 90 7  0°)^m  (to  scale) 
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Figure  4.  Buckling  Shapes  for  Case  1  with  p  -  0.0002646:  For  Case  1,  total 
volume  0  =  0.04,  with  ply  angle  orientation  of  (0 90°/  0°)I/B.  Re¬ 
presented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  5.  Geometry  for  Case  1  with  p  =  0.0002822:  For  Case  1,  total  volume  0 
=  0.04,  with  ply  angle  orientation  (0 °btaml  90°/  0°)^„  (to  scale) 

Of  the  two  optimum  solutions  found  for  Case  1,  the  global  optimum  occurs 
when  modes  1,  2,  and  4  are  simultaneous.  This  is  not  the  expected  result  in  which  all 
four  buckling  modes  occur  simultaneously.  The  increase  in  size  of  the  stiffener  for  the 
global  optimum  result  allows  the  t2  (0°)  layer  to  decrease  down  to  the  lower  limit.  This 
causes  the  stiffener  to  be  relatively  rig;d  and  the  plate  to  be  relatively  weak  in  the  loading 
direction.  Therefore,  for  the  global  optimum  case,  a  much  greater  portion  of  the  load 
is  carried  by  the  stiffener. 
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Figure  6.  Buckling  Shapes  for  Case  1  with  p  =  0.0002822:  For  Case  1,  total 
volume  0  =  0.04,  with  ply  angle  orientation  of  (0°^/  90°/  0°)vw.  Re¬ 
presented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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2.  Case  2:  0  =  0.04  and  Ply  Orientation  (0 °kmm!  0°/  90°)>x- 

For  this  configuration,  only  one  optimum  solution  is  identified.  It  is  listed  in 
Table  2  below.  The  same  results  are  presented  graphically  in  Figures  7  and  8.  This 
turns  out  to  be  the  best  possible  configuration  of  those  tested.  Note  that  here  all  four 
modes  buckle  simultaneously. 


Figure  7.  Geometry  for  Case  2  with  p-  -  0.0003949:  For  Case  2,  total  volume  0 


*  0.04,  with  ply  angle  orientation  (0°^/  0°/  90o)„„  (to  scale) 
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Figure  8.  Buckling  Shapes  for  Case  2  with  p  -  0.0003949:  For  Case  2,  total 
volume  0  *  0.04,  with  ply  angle  orientation  of  (0o*,w/  0°/  90%,,.  Re¬ 
presented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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3.  Case  3:  8  =  0.04  and  Ply  Orientation  (90 °tmml  90°/  0°),r. 

For  this  configuration,  two  different  optimum  solutions  are  identified.  They  are 
listed  in  Table  3  below.  The  same  results  are  presented  graphically  in  Figures  9-12. 


Table  3.  RESULTS:  8  =  0.04,  PLY  ANGLE  (90° 4£^/  90°/  0°)s 
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Figure  10.  Buckling  Shapes  for  Case  3  with  p  -  0.0000677:  For  Case  3,  total 
volume  ©  *  0.04,  with  ply  angle  orientation  of  (90° I  90°/  0°)IT„. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Of  the  two  optima  found,  the  global  optimum  occurs  when  three  buckling 
modes  occur  simultaneously.  No  optimum  was  located  that  gave  simultaneous  buckling 
with  all  four  modes.  For  this  case,  the  stifTener  fibers  are  in  the  90°  orientation  so  the 
stiffener  adds  much  less  strength  to  the  plate  than  when  the  stiffener  fibers  are  in  the  0° 
orientation.  Thus  more  of  the  total  volume  must  be  placed  into  the  r2  (0°)  lamina.  Be¬ 
cause  of  the  volume  of  material  in  the  t}  lamina  and  the  stifTener,  there  is  an  insufficient 
volume  of  material  to  allow  the  r,  (90°;  lamina  to  increase  from  the  lower  limit.  Because 
of  this,  it  is  impossible  to  obtain  simultaneous  buckling  with  all  four  modes. 

For  the  local  optima  case  where  p  =  0.0000677,  the  stiffened  plate  is  trying  to 
reduce  itself  down  to  a  single  lamina  with  0°  fiber  orientation.  From  this  configuration, 
any  incremental  increase  in  the  stiffener  size  causes  a  decrease  in  buckling  load.  By  ini- 
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tially  increasing  the  stiffener  size,  the  load  share  carried  by  the  stiffener,  pb,  increases  with 
the  increase  in  the  cross  sectional  area,  sy  The  initial  increase  in  pt  is  greater  than  the 
increase  in  stiffness  caused  by  the  increase  in  sb,  therefore,  the  buckling  load  initially 
decreases  creating  the  local  optimum. 


Figure  12.  Buckling  Shapes  for  Case  3  with  p  =  0.0001804:  For  Case  3,  total 
volume  0  =  0.04,  with  ply  angle  orientation  of  (90o*.m/  90°/  0°),r„. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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4.  Case  4:  0  =  0.04  and  Ply  Orientation  (90°^„/  0°/  90°)J>r- 

For  this  configuration,  three  different  optimum  solutions  are  identified.  They 
are  listed  in  Table  4  below.  The  same  results  are  presented  graphically  in  Figures  13-18. 


Table  4.  RESULTS;  0  =  0.04,  PLY  ANGLE  (90°,^/  0°/  90°)JK* 


h 

b 

h 

P 

Buckling 

Modes 

0.01900 

0.00100 

0.00100 

0.00100 

0.0000677 

1 

0.01636 

0.00100 

0.01845 

0.14210 

0.0001345 

1,2,4 

0.00778 

0.00816 

0.02830 

0.14310 

0.0001702 

1.2, 3.4 

Figure  13.  Geometry  for  Case  4  with  p  =  0.0000677:  For  Case  4,  total  volume 
©  =  0.04,  with  ply  angle  orientation  (90°*^*/  0°/  90o)^M  (to  scale) 
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Figure  14.  Buckling  Shapes  for  Case  4  with  p  =  0.0000677:  For  Case  4,  total 
volume  ©  =  0.04,  with  ply  angle  orientation  of  (90°*.„/  0°/  90°)1>m. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  15.  Geometry  for  Case  4  with  p  =  0.0001345:  For  Case  4,  total  volume 
0  =  0.04,  with  ply  angle  orientation  (90°^ /0°/  90°),y„  (to  scale) 

For  this  case,  the  global  optimum  occurs  when  all  four  buckling  modes  occur 
simultaneously.  However,  this  global  optimum  is  less  than  the  global  optimum  of  Case 
3.  It  appears  as  though  the  inner  lamina  (/2)  is  the  preferred  location  for  the  0°  fibers. 

The  configuration  which  gives  p  —  0.0001345  corresponds  closely  to  the  Case 
3  global  optimum.  The  only  difference  is  that  more  0°  fiber  orientation  volume  is  re¬ 
quired  when  it  is  located  at  the  outer  layer  (r,).  The  local  optimum,  p  =  0.0000677, 
corresponds  exactly  to  the  Case  3  local  optimum  with  the  same  load  value.  Here  the 
same  reasons  apply  for  the  existence  of  this  local  optimum. 


Figure  16.  Buckling  Shapes  for  Case  4  with  p  =  0.0001345:  For  Case  4,  total 
volume  0  =  0.04,  with  ply  angle  orientation  of  (90°*.„/  0°/  90°),^. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  17.  Geometry  for  Case  4  with  p  =  0.0001702:  For  Case  4,  total  volume 
0  =  0.04,  with  ply  angle  orientation  (90°*,,,/ 07  90°)^m  (to  scale) 


Figure  18.  Buckling  Shapes  for  Case  4  with  p  =  0.0001702:  For  Case  4,  total 
volume  ©  *  0.04,  with  ply  angle  orientation  of  (90°^/  0°/  90o)Iy„. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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5.  Case  5:  0  =  0.015  and  Ply  Orientation  (0 90°/  0°)„. 

For  this  configuration,  two  optimum  solutions  are  identified.  They  are  listed  in 
Table  5  below.  The  same  results  are  presented  graphically  in  Figures  19-22. 


Table  5.  RESULTS:  0  =  0.015,  PLY  ANGLE  (0°,^M/  90°/  0°)sy„ 


h 

b 

h 

P 

Buckling 

Modes 

0.00298 

0.00411 

0.01008 

0.04029 

0.0000156 

1,2.3, 4 

0.00576 

0.00100 

0.02720 

0.02720 

0.0000124 

1 

Figure  19.  Geometry  for  Case  5  with  p  =  0.0000156:  For  Case  5,  total  volume 
©  =  0.015,  with  ply  angle  orientation  (0°*„m/  90°/  0°)^M  (to  scale) 


Figure  20.  Buckling  Shapes  for  Case  5  with  p  =  0.0000156:  For  Case  5,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (0°*,„/ 90°/  0°),r„. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  21.  Geometry  for  Case  5  with  p  *  0.0000124:  For  Case  5,  total  volume 
0  =  0.015,  with  ply  angle  orientation  (0°*„„/  90°/  0°),r„  (to  scale) 


Here  the  global  optimum  is  found  when  the  four  buckling  modes  occur  simul¬ 
taneously.  This  global  optimum  corresponds  to  the  global  optimum  found  in  Case  1. 
This  follows  since  Case  5  is  the  same  as  Case  1,  except  for  the  smaller  volume  of  material 
available.  The  local  optimum  found  with  p  -  0.0000124  is  not  a  uniquely  determined 
solution.  The  total  stiffener  volume  is  uniquely  defined  but  not  the  individual  dimen¬ 
sions  of  the  stiffener  height  and  width.  This  occurs  because  the  first  buckling  mode  (the 
only  active  mode  in  this  case)  does  not  cause  any  flexure  of  the  stiffener.  The  active 
buckling  shape  only  twists  the  stiffener  so  that  only  its  volume  and  not  its  cross  sectional 
shape  is  important  to  the  overall  buckling  stiffness. 
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Figure  22.  Buckling  Shapes  for  Case  5  with  p  =  0.0000124:  For  Case  5,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (0°*^,/ 90°/ 0°)„m. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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6.  Case  6:  0  *  0.015  and  Ply  Orientation  (0 °Utml  0°/  90°),y. 

For  this  configuration,  only  one  optimum  solution  is  identified.  It  is  listed  in 
Table  6  below.  The  same  results  are  presented  graphically  in  Figures  23-24.  This  case 
is  similar  to  Case  2  except  here  there  is  insufficient  material  available  to  produce  simul¬ 
taneous  buckling  of  all  four  modes. 


Table  6.  RESULTS:  0  -  0.015,  PLY  ANGLE  (0°ftM,/  0°/  90°)jy* 


Figure  23.  Geometry  for  Case  6  with  p  =  0.0000179:  For  Case  6,  total  volume 


0  =  0.015,  with  ply  angle  orientation  (0°*,,,/  0°/  90°),,,,  (to  scale) 
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Figure  24.  Buckling  Shapes  for  Case  6  with  p  =  0.0000179:  For  Case  6,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (0°*„J  0°/  90°), 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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7.  Case  7:  0  =  0.015  and  Ply  Orientation  (90 90°/  0°)„. 

For  this  configuration,  two  different  optimum  solutions  are  identified.  They  are 
listed  in  Table  7  below.  The  same  results  are  presented  graphically  in  Figures  25-28. 


Table  7.  RESULTS:  0  =  0.015,  PLY  ANGLE  (90°,£,„/  90°/  0°)sm 


f, 

b 

h 

P 

Buckling 

Modes 

0.00100 

0.00650 

0.00100 

0.00100 

0.0000036 

1 

0.00100 

0.00584 

0.00669 

0.09761 

0.0000130 

1,2.4 
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Figure  26.  Buckling  Shapes  for  Case  7  with  p  —  0.0000036:  For  Case  7,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (90°^/ 90°/  0°),,.,. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  27.  Geometry  for  Case  7  with  p  -  0.0000130:  For  Case  7,  total  volume 
0  =  0.015,  with  ply  angle  orientation  (90°^/  90°/  0°)f>m  (to  scale) 


The  optima  of  Case  7  are  closely  related  to  the  optima  of  Case  3.  The  only 
difference  in  the  two  cases  is  the  total  volume  of  material  available.  Here  the  same  rel¬ 
ative  magnitude  between  the  two  different  optima  are  present  with  the  global  optima 
again  being  the  configuration  with  three  buckling  modes  occurring  simultaneously.  As 
in  Case  3,  here  too,  there  is  insufficient  material  available  to  obtain  the  case  where  all 
four  buckling  modes  occur  simultaneously. 
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Figure  28.  Buckling  Stupes  for  Case  7  with  p  -  0.0000130:  For  Case  7,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (90°*„„/ 90°/ 0°)„„. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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8.  Case  8:  0  =  0.015  and  Ply  Orientation  (90°^.  /  0°/  90°),,. 

For  this  configuration,  two  different  optimum  solutions  are  identified.  They  are 
listed  in  Table  8  below.  The  same  results  are  presented  graphically  in  Figures  29-32. 


Table  8.  RESULTS:  0  =  0.015,  PLY  ANGLE  (90 0°/90°)s 
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Figure  30.  Buckling  Shapes  for  Case  8  with  p  =  0.0000036:  For  Case  8,  total 
volume  ©  *  0.015,  with  ply  angle  orientation  of  0°/  90o),r„. 

Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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Figure  31.  Geometry  for  Case  8  with  p  -  0.0000125:  For  Case  8,  total  volume 
0  =  0.015,  with  ply  angle  orientation  (90°^/  0°/  90°)^m  (to  scale) 


As  in  Case  4,  the  global  optimum  solution  for  Case  8  occurs  when  all  four 
buckling  modes  occur  simultaneously.  Also,  there  is  a  similar  1  al  optimum  to  that  of 
Case  4,  which  has  only  one  active  buckling  mode.  This  case  is  present  for  the  same 
reasons  as  given  in  Case  3.  One  difference  between  Case  8  and  Case  4  is  that  no  local 
optimum  in  which  three  buckling  modes  occur  simultaneously  could  be  found  for  Case 
8  as  was  found  in  Case  4. 


48 


Figure  32.  Buckling  Shapes  for  Case  8  with  p  =  0.0000125:  For  Case  8,  total 
volume  0  =  0.015,  with  ply  angle  orientation  of  (90 °tnml  0°/  90°)rj^. 
Represented  above  are  the  buckling  shapes  of  the  simultaneous  overall 
plate  buckling  modes. 
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B.  OPTIMIZATION  EFFICIENCY 

To  compare  the  efficiencies  of  the  three  optimization  routines,  two  separate  criteria 
are  used.  These  criteria  are:  the  total  execution  time  required  to  converge  to  an  opti¬ 
mum  solution;  and  the  ability  of  an  optimizer  to  converge  to  a  solution  when  the  initial 
guess  is  far  from  the  solution  (i.e.  globally  convergent).  In  order  to  compare  the  relative 
efficiencies  of  the  three  optimization  methods,  the  exact  same  problem  is  solved  with 
each  of  the  optimization  schemes.  For  each  run,  the  same  initial  guess  is  used  and  the 
actual  execution  time  is  measured.  These  results  are  presented  in  Table  9  below. 


Table  9.  RUN  TIMES:  0  =  0.04,  PLY  ANGLE  (0%^*/  90°/  0°)SKM 


Opti¬ 

mization 

Program 

Initial  Guess 

Exe¬ 
cution 
Time  in 
CPU 
seconds 

h 

h 

b 

h 

P 

0.01000 

0.00100 

0.10000 

0.10000 

0.0001000 

Output 

h 

b 

h 

p 

IMSL 

0.01521 

0.00100 

0.05794 

0.06526 

0.0000282 

35.8 

DOT  1 

0.01526 

0.00100 

0.05831 

0.06510 

0.0000287 

72.9 

DOT  2 

0.01521 

0.00100 

0.05795 

0.06525 

0.0000282 

132.4 

From  these  results,  it  can  be  seen  that  for  this  specific  design  problem,  the  IMSL 
routine  is  clearly  the  most  efficient  in  terms  of  the  execution  time.  The  relative  speed 
of  convergence  for  the  IMSL  routine  is  due  to  the  high  degree  of  nonlinearity  in  the 
design  problem.  Because  the  IMSL  routine  uses  sequential  quadratic  programming  as 
its  search  algorithm,  which  is  based  on  a  second  order  derivative,  it  allows  each  iteration 
in  the  optimum  search  to  be  more  profitable  than  those  of  the  DOT  algorithms.  Both 
the  DOT  methods  (MFD  and  SLP)  have  search  directions  that  are  based  on  first  order 
derivatives.  Although  the  search  direction  calculation  for  both  DOT  methods  is  simpler, 
it  is  not  as  profitable  in  each  iteration  in  the  optimum  search  as  that  of  the  IMSL 
method.  Thus,  more  iterations  are  required  for  optimum  solution  convergence  with  the 
DOT  methods.  The  savings  made  in  the  calculation  time  for  each  search  direction  in  the 
DOT  code  is  not  nearly  enough  to  make  up  for  the  additional  calculation  time  required 
for  the  extra  iterations  in  the  optimization  Therefore,  due  to  the  highly  nonlin  na¬ 
ture  of  this  problem,  the  IMSL  routine  is  -  efficient  in  terms  of  execution  tin 
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The  second  measure  of  the  optimizer's  efficiency  is  its  ability  to  be  globally  conver¬ 
gent.  That  is,  to  converge  to  some  solution  from  almost  any  starting  point  [Ref.  32:  p. 
5J.  Since  the  location  and  number  of  optima  are  unknown  at  the  start  of  the  optirruza- 
tion  process,  the  initial  guesses  are  scattered  throughout  the  design  space  in  an  attempt 
to  identify  all  possible  optima.  Because  the  1MSL  routine  is  based  on  a  second  order 
search,  it  is  better  suited  to  converge  when  the  initial  guess  is  relatively  far  away  from 
an  optima.  Both  the  DOT  methods  were  unable  to  converge  to  an  optima  on  many  runs 
when  the  initial  guess  was  relatively  far  from  any  optima.  Qualitatively,  the  IMSL  rou¬ 
tine  produced  actual  optimum  convergence  ten  fold  more  often  than  either  of  the  DOT 
methods.  Therefore,  from  the  standpoint  of  globally  convergent  efficiency,  the  I  MSL 
routine  is  clearly  superior  for  this  specific  nonlinear  optimization  problem. 


V.  CONCLUSIONS 


Many  times  it  is  difficult  to  draw  general  conclusions  for  design  problems.  This  case 
is  no  exception,  however,  there  are  three  main  areas  where  specific  conclusions  for  this 
particular  design  problem  can  be  drawn.  These  areas  are:  the  general  character  of  the 
design  problem;  the  best  configuration  of  the  investigated  stacking  sequences;  and  the 
relative  efficiencies  of  the  three  optimization  routines  used. 

The  first  general  conclusion  drawn  is  the  character  of  the  design  problem  for  the 
maximum  buckling  load  of  the  stiffened  composite  plate.  The  design  optimization 
problem  turns  out  to  be  a  highly  nonlinear  problem.  For  almost  all  cases  investigated, 
there  were  multiple  optimum  solutions  indicating  a  nonlinear  solution  to  the  design  op¬ 
timization  problem.  Some  cases  had  as  many  as  three  different  optima.  Even  the  case 
where  only  one  optimum  solution  is  identified  cannot  be  assumed  to  be  an  indication 
of  a  linear  problem.  This  must  be  viewed  as  the  only  solution  identified  and  not  as  the 
only  solution  that  exists.  Additionally,  these  solutions  were  all  located  in  a  limited  de¬ 
sign  space.  One  would  expect  even  more  optima  to  be  found  by  increasing  the  size  of 
the  design  space.  Finally,  the  order  of  the  governing  partial  differential  equation  would 
lead  one  to  assume  a  highly  nonlinear  problem  without  examining  the  results.  The  re¬ 
sults  merely  confirm  this  assumption. 

The  second  area  of  general  conclusions  is  that  of  the  best  design  configuration.  For 
both  total  volume  constraints  (0  =  0.04  and  0.015)  examined,  the  best  stacking  config¬ 
uration  was  the  (0°4mw/  0°/  90°)„m  used  in  cases  2  and  6.  The  use  of  the  0°4„„  stiffener 
fiber  orientation  was  always  superior  to  the  90°^  stiffener  fiber  orientation  regardless 
of  the  laminae  ply  orientation.  This  led  to  the  conclusion  that  the  0 aknm  stiffener  fiber 
orientation  is  the  most  important  factor  in  the  design  optimization.  The  laminae  fiber 
orientation  of  (0°/90°)OM  was  best  when  used  with  the  90°*,m  stiffener  fiber  orientation. 
However,  when  used  with  the  0°^  stiffener  fiber  orientation,  the  {90°  IQ°)^m  laminae  fi¬ 
ber  orientation  was  best.  This  indicates  that  for  best  results,  the  inner  lamina  ( i }  layer) 
should  have  its  fiber  orientation  perpendicular  to  that  of  the  stiffener  fiber  orientation. 

The  third  area  of  general  conclusions  is  the  relative  efficiencies  of  the  three  opti¬ 
mization  routines  used.  Because  of  the  highly  nonlinear  nature  of  this  specific  design 
optimization  problem,  the  IMSL  subroutine  DNCONG  proved  to  be  the  most  efficient, 
both  in  terms  of  convergence  speed  and  in  its  ability  to  be  globally  convergent  (i.e. 
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converge  to  an  optimum  from  almost  any  starting  point).  Because  the  IMSL  routine  is 
based  on  a  second  order  search  direction,  it  is  better  suited  to  optimize  the  highly  non¬ 
linear  design  problem.  In  contrast  to  this,  both  of  the  DOT  methods  have  search  di¬ 
rections  that  are  based  on  first  order  search  directions  that  make  them  ill-suited  for  such 
a  highly  nonlinear  optimization  problem. 

It  is  important  to  note  that  the  conclusions  of  this  research  are  specific  to  the  spe¬ 
cific  design  optimization  problem  investigated.  Generalizations  of  these  results  to  other 
problems  must  be  examined  closely  due  to  the  nonlinear  nature  of  the  problem.  Any 
design  optimization  problem  of  a  nonlinear  system  must  be  closely  scrutinized  by  the 
designer  to  ensure  that  all  important  aspects  of  the  design  are  considered.  Computer- 
aided  design  optimization  can  be  a  very  powerful  tool  for  a  designer,  however,  caution 
must  be  exercised  to  ensure  that  a  local  optimum  design  is  not  selected  as  the  best  pos¬ 
sible  design  without  investigating  the  entire  design  space  for  other  optima. 
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APPENDIX  A.  DESIGN  ANALYSIS  SUBROUTINES 


SUBROUTINE  EIGANL4 


LINEAR  FINITE  ELEMENT  ANALYSIS  OF  LAMINATED  COMPOSITE  PLATE 
(16  DOF  RECTANGULAR  PLATE  ELEMENT) 

3  LOWEST  EIGENVLAUES  AND  THEIR  DERIVATIVES  WRT  DESIGN  VARIABLES 


IMPLICIT  REAL*8(A-H,0*Z) 

COMMON / ELEMNT /  EK(16,16,6),  EKT(16,16),  EG(16,16) 

COMMON/ELEMNB/  ARSTF(4,4),  ARGEM(4,4) 

COMMON /PROPT1/  El(15),  E2(15),  G12(15),  P012(15),  ANGL( 15) ,  Z( 16) 

*  ,  NOD( 150 ,4) ,  NX,  NY,  NSDF,  LSDF(150) 

COMMON/PROPT2/  DM(3,3) 

COMMON / PROPT3 /  EBEAM,  ABEAM,  BMINER ,  EPLT,  APLT ,  THICK,  NODB(18,2) 
COMMON / PR0PT4/  HT(15),  AAA,  BBB ,  NLYR ,  NLYR1,  NLYR2 ,  NLYR3 
COMMON/TOTAL  /  GK( 150 , 150) ,GG( 150 , 150) 

COMMON/BAND/  BANDK(45 , 150) ,  BANDKR(45 , 150) ,  BANDG(45 , 150) 
COMMON/EIG/NBAND ,  NREQ,  NEQ 

C0MM0N/EIG2/  EGVC( 150,5),  EGVCO( 150 ,5) ,  EGVL(5,4),  EGVL0(5,4), 

*  DTHK(50 ,5) ,  DGD(5),  EGVCB( 150 ,5) ,  IDMOD 
DIMENSION  DBANDK(45 , 150) ,  WK(150),  WK1(150),  WK2(150) 

COMMON/BEAM/  E1B,  E2B,  G12B,  P012B,  BBM,  HBM 
COMMON/ BEAM2/  DBTHK(50) ,  PBEAM 
COMMON/NORMAL/  VNORMAL 

CALL  PROPTY 
CALL  ELESTF 
CALL  ASSEMB 
CALL  EIGEN 
IDMOD  =  1 

DO  11  1=1,3 

EGVL(I,1)  =  EGVL( 1,1) /VNORMAL 
11  CONTINUE 

---  SAVE  K,  G 

NBAND1  =  NBAND  +  1 
DO  150  1=1,  NREQ 

DO  140  J=l,  NBAND 1 

DBANDK(J.I)  =  BANDKR(J.I)  -  EGVL( 1 , 1)*BANDG( J , I ) 

140  CONTINUE 

150  CONTINUE 
C 
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c  === 
c 


1305 


1320 


1307 


2320 

C 

C  *** 
C 


4320 

C 

C 

C 


C 

c  --- 
c 


280 

300 

C 


OBTAIN  DERIVATIVES  OF  P  WRT  T 

DO  1305  J=l,  NREQ 
WK( J)  =  EGVC(J, 1) 

CONTINUE 
LC  =  45 

CALL  BNDMUL( BANDG , UK , NBAND 1 , NREQ , WK2 , LC ) 

DGD(l)  =  0 . DO 
DO  1320  J=l,  NREQ 

DGD(l)  =  DGD(l)  +  WK(J)*WK2(J) 

CONTINUE 

IF(IDMOD.GE.l)  THEN 
DO  1307  J=l,  NREQ 
WK1(J)  =  EGVC( J ,2) 

CONTINUE 
LC  =  45 

CALL  BNDMUL( BANDG , WK 1 , NBAND 1 , NREQ , WK2 , LC ) 

DGD(2)  =  0 . DO 
DO  2320  J=l,  NREQ 

DGD(2)  =  DGD(2)  +  WK1(J)*WK2( J) 

CONTINUE 

FOR  3RD  EIGENVECTOR 

CALL  BNDMUL ( BANDG , EGVC ( 1 , 3 ) , NB AND 1 , NREQ , WK2 , LC ) 
DGD( 3)  =  0 . DO 
DO  4320  J=l,  NREQ 

DGD( 3)  =  DGD( 3)  +  EGVC( J  ,  3)*WK2( J) 

CONTINUE 
END  IF 


DEL  =  1.0D-05 
DO  380  IV=1,  NLYR2 

IF( IV . EQ . NLYR1)  THEN 
BBM  =  BBM  +  DEL 
ELSE  IFCIV.EQ.NLYR2)  THEN 
HBM  =  HBM  +  DEL 
ELSE 

HT(IV)  =  HT(IV)  +  DEL 
END  IF 
CALL  PROPTY 
CALL  ELESTF 
CALL  ASSEMB 

GET  DK 

DO  300  J=l,  NREQ 

DO  280  K=1 ,  NBAND1 

BANDKR(K, J)  =  (BANDKR(K.J)  -  EGVL( 1 , 1)*BANDG(K,J) 
+  -  DBANDK(K,J))/DEL 

CONTINUE 
CONTINUE 
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c  -- 
c 


320 

C 

C  --- 
C 


3320 

C 


5320 


C 

c  --- 
c 


380 

C 


4694 

C 


1111 

C 


C 

C 

C 

C 

C 


GET  U*DK*U 
LC  =  45 

CALL  BNDMUL( BANDKR , WK , NBAND1 , NREQ ,  WK2 , LC) 
DTHKCIV, 1)  =  0 .DO 
DO  320  J=l,  NREQ 

DTHKCIV, 1)  =  DTHKCIV, 1)  +  WK(J)*WK2CJ) 
CONTINUE 

DTHKC IV ,  1)  =  DTHK(IV,1)/DGD(1) 


GET  U*DK*U 


IF(IDMOD.GE.l)  THEN 
LC  =  45 

CALL  BNDMULC  BANDKR , WK 1 , NBAND 1 , NREQ , WK2 , LC ) 

DTHK( IV ,2)  =  0 . DO 
DO  3320  J=l,  NREQ 

DTHK( IV, 2)  =  DTHK( IV, 2)  +  WK1( J)*WK2( J) 
CONTINUE 

DTHK( IV ,2)  =  DTHKCIV ,2)/DGD(2) 

CALL  BNDMULC BANDKR ,EGVC( 1 ,3) .NBANDl ,NREQ ,WK2 ,LC) 
DTHKCIV, 3)  =  0.D0 
DO  5320  J=l,  NREQ 

DTHKCIV, 3)  =  DTHKCIV, 3)  +  EGVCCJ,3)*WK2CJ) 
CONTINUE 

DTHKCIV,))  =  DTHKC IV , 3)/DGDC  3) 

END  IF 

ORIGINAL  HT 

IF(IV.EQ.NLYRl)  THEN 
BBM  =  BBM  -  DEL 
ELSE  IFCIV.EQ.NLYR2)  THEN 
HBM  =  HBM  -  DEL 
ELSE 

HT(IV)  =  HTCIV)  -  DEL 
END  IF 
CONTINUE 

DO  4694  I=1,NLYR2 
DO  4694  J=1 , 3 

DTHKC I, J)  =  VNORMAL*DTHK(I,J) 

CONTINUE 

DO  1111  1=1,3 

EGVLCI.l)  =  EGVLCI,l)*VNORMAL 
CONTINUE 

RETURN 

END 
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SUBROUTINE  INPUTB 


C 
C 

C 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/PROPTO/  THETA,  ALEN ,  BLEN 

COMMON/PROPT1/  El(15),  E2(15),  G12(15),  P012(15),  ANGL(IS),  Z( 16) 

2  ,NOD( 150,4),  NX,  NY,  NSDF ,  LSDF(150) 

C0MM0N/PR0PT4/  HT(15),  AAA,  BBB ,  NLYR ,  NLYR1 ,  NLYR2 ,  NLYR3 
COMMON/BEAM/  E1B ,  E2B ,  G12B,  P012B,  BBM ,  HBM 
C 

C  READ  IN  GEOMETRIC  PARAMETERS  AND  PHYSICAL  PROPERTIES 
C 

READ(5 ,*)  ALEN,  BLEN 
READ(5 ,*)  NLYR, THETA 
NLYR1  =  NLYR  +  1 

NLYR 2  =  NLYR  +  2 

NLYR3  =  NLYR  +  3 

DO  150  1=1,  NLYR 

READ(5 ,*)  E1(I),  E2( I )  ,  G12(I),  P012(I),  ANGL(I) 

150  CONTINUE 

READ( 5 ,*)  NX,  NY 
READ( 5 ,*)  E1B,  E2B ,  G12B,  P012B 
RETURN 
C 

C  _ 

c 

SUBROUTINE  INPUTC 

C  _ 

c 

IMPLICIT  REAL*8(A-H ,0-Z) 

EXTERNAL  GRMESH 

COMMON/PROPTO/  THETA,  ALEN,  BLEN 

COMMON/PROPT1/  El(15),  E2(15),  G12(15),  P012(15),  ANGL( 15) ,  Z( 16) 

*  ,  NOD( 150 ,4) ,  NX,  NY,  NSDF,  LSDF(150) 

COMMON/PROPT3/  EBEAM,  ABEAM,  BMINER,  EPLT,  APLT,  THICK,  NODB(18,2) 
C0MM0N/PR0PT4/  HT(15),  AAA,  BBB,  NLYR,  NLYR1 ,  NLYR2 ,  NLYR3 
COMMON/EIG/NBAND ,  NREQ ,  NEQ 

COMMON/EIG2/  EGVC(150,5),  EGVCOC 150,5) ,  EGVL(5,4),  EGVLO(5,4), 

*  DTHK(50 , 5) ,  DGD(5) ,  EGVCBC 150,5) ,  IDMOD 
COMMON/BEAM/  E1B ,  E2B,  G12B,  P012B,  BBM,  HBM 

C 

Z(l)  =  0 . DO 
DO  150  1=1,  NLYR 

Z(I+1)  =  Z(I)  +  HT( I ) 

150  CONTINUE 

THICK  =  Z(NLYR1)*2 . 0D0 
APLT  =  BLEN*THICK 
C 

NEL  =  NX*NY 

CALL  GRMESH(NNM ,NEM) 

C  ***  BEAM  CONNECTIVITY  MATRIX  *** 

C 

NX1  =  NX  +  1 
NY1  =  NY  +  1 
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124 

C 

c  *** 
c 


200 


250 


270 


300 

400 

C 


C 

C  -- 
C 


721 

722 

C 

C 

C 

C 


NODO  =  NY/2*NX1  +  1 
DO  124  1=1,  NX 
NODN  =  NODO  +  1 
NODB(I.l)  =  NODO 
NODB( 1,2)  =  NODN 
NODO  =  NODN 
CONTINUE 

BOUNDARY  CONDITIONS  OF  SIMPLE  SUPPORT  *** 

NSDF  =  12  +  4*(NX-1)  +  4*(NY-1) 

N  =  0 

DO  400  K=1 ,  NY1 
DO  300  J=l,  NX1 

IF(K.NE. 1.AND.K.NE.NY1)  THEN 

IF(J.NE.l.AND.J.NE.NXl)  THEN 
GO  TO  300 

ELSE 

DO  200  1=1,  2 
N  =  N  +  1 

LSDF(N)  =  4*((K-1)*NX1  +  J  -  1)  +  I 
CONTINUE 
END  IF 

ELSE 

IF(J.NE.l.AND.J.NE.NXl)  THEN 
DO  250  1=1,  3,  2 
N  =  N  +  1 

LSDF(N)  =  4*((K-1)*NX1  +  J  -  1)  +  I 
CONTINUE 

ELSE 

DO  270  1=1,  3 
N  =  N  +  1 

LSDF(N)  =  4*((K-1)*NX1  +  J  -  1)  +  I 
CONTINUE 
END  IF 
END  IF 
CONTINUE 
CONT .NUE 


AAA  =  ALEN/NX 

BBB  =  BLEN/NY 

NEQ  =  4*(NX+1)*(NY+1) 

NREQ  =  NEQ  -  NSDF 

CLEAR  EGVCB 

DO  722  1=1,  3 

DO  721  J=l,  NREQ 
EGVCB(J.I)  =  0.D0 
CONTINUE 
CONTINUE 
RETURN 
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SUBROUTINE  PROPTY 


C 

C 


IMPLICIT  REAL*8(A-H ,0-Z) 

COMMON/ PR0PT4/  HT(15),  AAA,  BBB ,  NLYR ,  NLYR1 ,  NLYR2 ,  NLYR3 
COMMON / PROPT 1 /  E1(1S),  E2(15),  G12(15),  P012(15),  ANGL( 15) ,  Z( 16) 

*  ,  N0D( 150,4),  NX,  NY,  NSDF,  LSDF(150) 

C0MM0N/PR0PT2/  DM(3,3) 

COMMON/ PR0PT3/  EBEAM,  ABEAM,  BMINER ,  EPLT,  APLT,  THICK,  NODB(18,2) 
COMMON/BEAM/  E1B ,  E2B,  G12B,  P012B,  BBM ,  HBM 
DIMENSION  DANGL( 15) 

Z(l)  =  0 . DO 
DO  5  1=1,  NLYR 

Z(I+1)  =  Z( I)  +  HT( I) 

5  CONTINUE 

DO  10  1=1,  3 
DO  10  J=l,  3 
DM( I , J)=0 . 0 
10  CONTINUE 

AM 11  =  0 . DO 
AM12  =  0 . DO 
AM13  =  0 . DO 
AM22  =  0.D0 
AM23  =  0.D0 
AM33  =  0.D0 

***  DO-LOOP  FOR  NUMBER  OF  LAYERS 
DO  30  1=1, NLYR 

***  CALCULATION  OF  THE  REDUCED  STIFFNESSES 

P021=E2( I )*P012( I ) /El( I ) 

BPO=l . 0-P012( I )*P021 
Q1 1=E1( I )/BP0 
Q12=P012( I)*E2( I)/BPO 
Q22=E2( I )/BPO 
Q66=G12( I) 

***  CALCULATION  OF  THE  TRANSFORMED  REDUCED  STIFFNESSES 

DANGL( I )=3 . 141592D0*ANGL( I )/ 180 . DO 
SC=DSIN(DANGL( I ) )*DCOS(DANGL( I ) ) 

S2=(DSIN(DANGL( I) ) )**2 
C2=(DCOS(DANGL( I) ) )**2 

QB 1 l=Qll*C2*C2+2 . 0*(Q12+2 . 0*Q66)*SC*SC+Q22*S2*S2 
QB22=Q11*S2*S2+2.0*(Q12+2.0*Q66)*SC*SC+Q22*C2*C2 
QB12=(Q11+Q22-4.0*Q66)*SC*SC+Q12*(S2*S2+C2*C2) 
QB16=(Q11-Q12-2.0*Q66)*SC*C2+(Q12-Q22+2.0*Q66)*SC*S2 
QB26=(Q1 1-Q12-2 . 0*Q66)*SC*S2+(Q12-Q22+2 . 0*Q66)*SC*C2 
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QB66=(Q11+Q22-2.0*Q12-2.0*Q66)*SC*SC+Q66*(S2*S2+C2*C2) 


***  CALCULATION  OF  THE  EXTENSIONAL  STIFFNESSES 

AMI 1=QB1 1*( (Z( 1+1) )-(Z( I ) ) )*2 . 0D0  +  AM11 
AM12=QB12*((Z(I+1))-(Z(I)))*2.0D0  +  AM12 
AM22=QB22*((Z(I+1))-(Z(I)))*2.0D0  +  AM22 
AM13=QB16*((Z(I+1))-(Z(I)))*2.0D0  +  AM13 
AM33=QB66*((Z(I+1))-(Z(I)))*2.0DO  +  AM33 
AM23=QB26*( (Z( I+1))-(Z(I)))*2. 0D0  +  AM23 

***  CALCULATION  OF  THE  BENDING  STIFFNESSES 

DM(1,1)=QB11*((Z(I+1))**3-(Z(I))**3)/1.5D0  +  DM(1,1) 
DM(1,2)=QB12*((Z(I+1))**3-(Z(I))**3)/1.5D0  +  DM(1,2) 
DM(2 , 2)=QB22*( (Z( 1+1 ) )**3- (Z( I ) )**3)/ 1 . SDO  +  DM(2,2) 
DM(1,3)=QB16*((Z(I+1))**3-(Z(I))**3)/1.5D0  +  DM(1,3) 
DM( 3 , 3)=QB66*( (Z( 1+1) )**3- (Z( I ) )**3)/ 1 . SDO  +  DM(3,3) 
DM(2,3)=QB26*((Z(I+1))**3-(Z(I))**3)/1.SD0  +  DM(2,3) 
30  CONTINUE 
AM21  =  AM 12 
AM31  =  AM13 
AM32  =  AM23 

***  FIND  E  OF  THE  PLATE  *** 

BLEN  =  BBB*NY 

THICK  =  Z(NLYR1)*2 . 0D0 

APLT  =  BLEN*THICK 


DENO  =  AM 1 1*AM22*AM33  +  AM 12*AM23*AM3 1*2 . DO  -  AM13**2*AM22  - 
*  AM12**2*AM33  -  AM11*AM23**2 

EPLT  =  DENO/ (AM22*AM33  -  AM23**2)/THICK 
DM( 2,1)  =  DM( 1,2) 

DM( 3,1)  =  DM( 1,3) 

DM( 3,2)  =  DM( 2 , 3) 

ABEAM  =  2 . D0*BBM*HBM 

EBEAM  =  E1B 

THIK2  =  THICK/2. DO 

BMINER  =  2 . D0*B8M/3 . D0*( (THIK2+HBM)**3-THIK2**3) 

RETURN 

END 


SUBROUTINE  ELESTF 


IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  GAUSS4(4) ,  GAUSS3(.3),  WT4(4),  WT3(3),  BG(16) 
COMMON/ELEMNT/  EK(16,16,6),  EKT(16,16),  EG(16,16) 

COMMON/ELEMNB/  ARSTF(4,4),  ARGEM(4 ,4) 

COMMON/PROPT2/  DM(3,3) 

COMMON/PROPT3/  EBEAM,  ABEAM,  BMINER,  EPLT,  APLT,  TT  'CK,  NODB(18,2) 
C0MM0N/PR0PT4/  HT(15),  AAA,  BBB ,  NLYR,  NLYR1,  NLYF  NLYR3 
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C 


DIMENSION  BK( 16,3) ,  WK1616( 16, 16) ,  WK163(16,3),  MK316(3,16), 
*  WK33(3,3) ,  DV(6) 

DATA  GAUSS4/ . 06943I8D0 ,0 . 3300095D0 ,0 . 6699905D0 ,0 . 9305681D0/ 
DATA  GAUSS3/ . 1 127017D0 , 0 . 5D0 , 0 . 88729833D0/ 

DATA  WT4/0 . 1739274D0 , 0 . 3260725D0 , 0 . 3260725D0 , 0 . 1739274D0/ 
DATA  WT3/0 . 2777777D0 ,0 . 4444444D0 , 0 . 2777777D0/ 

***  INTERPOLATION  FUNCTIONS  AND  THEIR  DERIVATIVES  *** 

FN1(P)  =  1 .DO  -  3.D0*P*P  +  2.D0*P*P*P 
FN2( P)  =  ( 1 . DO  -  2.D0*P  +  P*P)*P 
FN3( P)  =  (3. DO  -  2.D0*P)*P*P 
FN4(P)  =  P*P*(1.D0  -  P) 

DFN1(P)  *  P*( -6 . DO  +  6.D0*P) 

DFN2(P)  =  1 . DO  -  4.D0*P  +  3.D0*P*P 
DFN3(P)  =  6.D0*P*(1.D0  -  P) 

DFN4(P)  =  P*(2 .DO  -  3 . D0*P) 

DDFNl(P)  =  -6. DO  +  12.D0*P 
DDFN2( P)  =  -4. DO  +  6.D0*P 
DDFN3( P)  =  6. DO  -  12.D0*P 
DDFN4(P)  =  2. DO  -  6.D0*P 

***  COMPONENT  OF  D  MATRIX 

DV( 1 )  =  DM(1,1) 

DV( 2)  =  DM( 1,2) 

DV(  3)  =  DM(  1,3) 

DV(4)  =  DM(2,2) 

DV(5 )  =  DM( 2,3) 

DV(6)  =  DM(3 , 3) 


***  INITIALIZE  THE  ELEMENT  MATRICES  AND  VECTORS 

DO  20  J=l,  16 
DO  20  1=1,  16 
EG( I , J)  =  0 .DO 
EKT( I , J)  =  0.D0 
20  CONTINUE 

DO  30  K=1 ,  6 

DO  30  J=l,  16 
DO  30  1=1,  16 

EK( I , J,K)  =  0 . DO 

30  CONTINUE 

***  NUMERICAL  INTEGRATION  BY  GAUSS  QUADRATURE 

DO  500  NI=1 , 3 
XI=GAUSS3(NI ) 

DXF1  =  DFNl(XI) 

DXF2  =  DFN2(XI ) 
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DXF3  =  DFN3(XI) 
DXF4  =  DFN4(XI ) 

DO  500  NJ=1 ,4 
ETA=GAUSS4(NJ) 
EF1  =  FNl(ETA) 
EF2  =  FN2(ETA) 
EF3  =  FN3(ETA) 
EF4  =  FN4(ETA) 


C 

BG( 1)  =  DXF1*EF1/AAA 
BG(2)  =  DXF1*EF2*BBB/AAA 
BG(3)  =  -DXF2*EF1 
BG(4)  =  DXF2*EF2*BBB 
C 

BG(5)  =  DXF3*EF1/AAA 
BG( 6)  =  DXF3*EF2*BBB/AAA 
BG( 7 )  =  DXF4*EF1 
BG(8)  =  -DXF4*EF2*BBB 
C 

BG(9)  =  DXF3*EF3/AAA 
BG(IO)  =  -DXF3*EF4*BBB/AAA 
BG(ll)  =  DXF4*EF3 
BG( 12)  =  DXF4*EF4*BBB 
C 

BG( 13)  =  DXF1*EF3/AAA 
BG( 14)  =  -DXF1*EF4*BBB/ AAA 
BG( 15)  =  -DXF2*EF3 
BG( 16)  =  -DXF2*EF4*BBB 
C 

DO  250  J=l,  16 
DO  250  1=1,  16 

EG( I , J)  =  EG(I.J)  +  WT3(NI)*WT4(NJ)*BG( I )*BG( J) 
250  CONTINUE 

C 

500  CONTINUE 


***  MULTIPLY  THE  LOAD  FACTOR  *** 


FACTOR  =  EPLT*THICK/(EPLT*APLT  +  EBEAM*ABEAM)*AAA*BBB 
DO  555  J=l,  16 
DO  555  1=1,  16 

EG( I , J)  =  EG( I , J)*FACTOR 
CONTINUE 


DO  600  NI=1 ,  4 
XI=GAUSS4(NI) 

XF1  =  FNl(XI) 

XF2  =  FN2(XI) 

XF3  =  FN3(XI) 

XF4  =  FN4(XI) 

DXF1  =  DFNl(XI) 
DXF2  =  DFN2(XI) 
DXF3  =  DFN3(XI ) 
DXF4  =  DFN4(XI ) 
DDXF1  =  DDFNl(XI) 
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DDXF2  =  DDFN2(XI) 

DDXF3  =  DDFN3(XI) 

DDXF4  =  DDFN4(XI) 

DO  600  NJ=1,  4 
ETA=GAUSS4(NJ) 

EF1  =  FNl(ETA) 

EF2  =  FN2(ETA) 

EF3  =  FN3(ETA) 

EF4  =  FN4(ETA) 

DEF1  =  DFNl(ETA) 

DEF2  =  DFN2(ETA) 

DEF3  =  DFN3(ETA) 

DEF4  =  DFN4(ETA) 

DDEF1  =  DDFNl(ETA) 

DDEF2  =  DDFN2(ETA) 

DDEF3  =  DDFN3(ETA) 

DDEF4  =  DDFN4(ETA) 

C 

BK( 1,1)  =  DDXF 1*EF 1 /AAA/ AAA 
BK(2 , 1)  =  DDXF1*EF2*BBB/AAA/AAA 
BKC3.1)  =  -DDXF2*EF1/AAA 
BK(4, 1)  =  DDXF2*EF2*BBB/AAA 
C 

BK( 5,1)  =  DDXF3*EF 1/ AAA/AAA 
BK( 6 , 1)  =  DDXF3*EF2*BBB/ AAA/AAA 
BK( 7 , 1)  =  DDXF4*EF1/AAA 
BK(8, 1)  =  -DDXF4*EF2*BBB/AAA 
C 

BK(9 , 1)  =  DDXF3*EF3/AAA/AAA 
BK( 10,1)  =  -DDXF3*EF4*BBB/AAA/AAA 
BK( 11,1}  =  DDXF4*EF3/AAA 
BK( 12,1)  =  DDXF4*EF4*BBB/AAA 
C 

BK( 13,1)  =  DDXF 1*EF3/ AAA/ AAA 
BK( 14,1)  =  -DDXF 1*EF4*BBB/ AAA/ AAA 
BK( 15,1)  =  -DDXF2*EF3/AAA 
BK( 16,1)  =  -DDXF2*EF4*BBB/AAA 
C 

BK(  1,2)  =  XF1*DDEF1/BBB/BBB 
BK(  2,2)  =  XF1*DDEF2/BBB 
BK(  3,2)  =  -XF2*DDEF1*AAA/BBB/BBB 
BK(  4,2)  =  XF2*DDEF2*AAA/BBB 
C 

BK(  5,2)  =  XF3*DDEF1/BBB/BBB 
BK(  6,2)  =  XF3*DDEF2/BBB 
BK(  7,2)  =  XF4*DDEF1*AAA/BBB/BBB 
BKC  8,2)  =  -XF4*DDEF2*AAA/BBB 
C 

BK(  9,2)  =  XF3*DDEF3/BBB/BBB 
BK( 10,2)  =  -XF3*DDEF4/BBB 
BK( 11,2)  =  XF4*DDEF3*AAA/BBB/BBB 
BK( 12,2)  =  XF4*DDEF4*AAA/BBB 
C 

BK( 13,2)  =  XF1*DDEF3/BBB/BBB 

BK( 14 , 2)  =  -XF1*DDEF4/BBB 

BK( 15,2)  =  -XF2*DDEF3*AAA/BBB/BBB 
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BK ( 16,2)  =  - XF 2*DDEF4*AAA / B  BB 
C 

BK(  1,3)  =  2 . DO*DXF 1*DEF 1/ AAA/ BBB 

BK(  2,3)  =  2 . D0*DXF 1*DEF2/ AAA 

BK(  3,3)  =  -2.D0*DXF2*DEF1/BBB 
BK(  4,3)  =  2 . D0*DXF2*DEF2 

C 

BK(  5,3)  =  2 .D0*DXF3*DEF1/AAA/BBB 

BK(  6,3)  =  2.D0*DXF3*DEF2/AAA 
BK(  7,3)  =  2 .D0*DXF4*DEF1/BBB 

BK(  8,3)  =  -2 . D0*DXF4*D£F2 
C 

BK(  9,3)  =  2 .D0*DXF3*DEF3/AAA/BBB 

BK( 10,3)  =  -2.D0*DXF3*DEF4/AAA 
BK( 11 , 3)  =  2.D0*DXF4*DEF3/BBB 
BK( 12,3)  =  2 . D0*DXF4*DEF4 

C 

BK( 13 ,3)  =  2.D0*DXF1*DEF3/AAA/BBB 
BK(14,3)  =  -2.D0*DXF1*DEF4/AAA 
BK( 15 ,3)  =  -2.D0*DXF2*DEF3/BBB 
BK( 16,3)  =  -2 . D0*DXF2*DEF4 
C 

C  ***  GET  TRANSPOSE  OF  BK 

C 

DO  570  1=1,  3 
DO  570  J=l,  16 

WK316( I , J)  =  BK( J , I ) 

570  CONTINUE 

C 

DO  580  KK=1 ,  6 
DO  560  JL*1,  3 
DO  560  IL=1 ,  3 

WK33(IL,JL)  =  0 . DO 
560  CONTINUE 

IF(KK.EQ.l)  THEN 
WK33( 1,1)  =  1 . DO 
ELSE  IF(KK.EQ.2)  THEN 
WK33( 1,2)  =  1 . DO 
WK33( 2 , 1)  =  1 .DO 
ELSE  IF(KK . EQ . 3)  THEN 
WK33( 1,3)  =  1 .DO 
WK33(3, 1)  =  1 . DO 
ELSE  IF(KK . EQ . 4)  THEN 
WK33( 2,2)  =  1 . DO 
ELSE  IF(KK . EQ . 5)  THEN 
WK33( 2 ,3)  =  1 .DO 
WK33( 3,2)  =  1 . DO 
ELSE  IF(KK .EQ . 6)  THEN 
WK33(3,3)  =  1 . DO 
END  IF 

CALL  MULTP2C BK , WK33 ,16,3, 3 , WK163) 

CALL  MULTP2(WK163 ,WK316 , 16,3,16,WK1616) 

C 

C  ***  MULTIPLY  WIGHT  FC .  AND  SUM  UP  FOR  ELEMENT  STIFF.  MATRIX 

C 

DO  550  J=l,  16 
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DO  550  1=1,  16 

EK(I,J,KK)  =  EK( I , J,KK)  +  WT4(NI)*WT4(NJ)*WK1616(I , J) 
550  CONTINUE 

580  CONTINUE 

C 

600  CONTINUE 
C 

C  ***  ASSEMBLE  FOR  ELEMENT  STIFFNESS  MATRIX 
C 

DO  620  1=1,  6 
DO  615  J=l,  16 
DO  615  K=1 ,  16 

EK(K , J , I)  =  EK(K,J,I)*AAA*BBB 
615  CONTINUE 

DO  620  J=l,  16 
DO  620  K=1 ,  16 

EKT(K.J)  =  EKT(K.J)  +  DV(I)*EK(K,J,I) 

620  CONTINUE 
C 

El  =  EBEAM*BMINER 
ARSTF(l.l)  =  EI*12.D0/AAA**3 
ARSTF( J ,2)  =  -EI*6.D0/AAA**2 
ARSTF( 1,3)  =  -ARSTFC 1,1) 

ARSTF( 1 ,4)  =  ARSTF( 1,2) 

ARSTF(2,2)  =  EI*4. DO/AAA 
ARSTFC2 ,3)  =  -ARSTF( 1,2) 

ARSTF(2 ,4)  =  EI*2. DO/AAA 
ARSTF( 3 ,3)  =  ARSTF(l.l) 

ARSTF(3,4)  =  -ARSTFC 1,2) 

ARSTF(4 ,4)  =  ARSTF( 2 ,2) 

C 

ARSTF(2 , 1)  =  ARSTF( 1,2) 

ARSTF(3, I)  =  ARSTFC 1,3) 

ARSTFC 3, 2)  =  ARSTF(2,3) 

ARSTFC 4 , 1)  =  ARSTFC 1,4) 

ARSTFC 4, 2)  =  ARSTFC 2, 4) 

ARSTFC4.3)  =  ARSTFC 3, 4) 

C 

C 

FACTOR  =  EBEAM*ABEAM/CEPLT*APLT  +  EBEAM*ABEAM) 

ARGEMC 1,1)  =  FACTOR* 1.2D0/ AAA 
ARGEMC 1,2)  =  -FACTOR*0 . IDO 
ARGEMC 1,3)  =  -ARGEMC 1,1) 

ARGEMC 1,4)  =  ARGEMC 1,2) 

ARGEMC2.2)  =  FACTOR*2 . DO/ 15 . D0*AAA 
ARGEMC 2, 3)  =  -ARGEMC 1,2) 

ARGEMC 2, 4)  =  -FACTOR*AAA/30 . DO 
ARGEMC 3, 3)  =  ARGEMC 1,1) 

ARGEMC3.4)  =  -ARGEMC 1,2) 

ARGEMC4.4)  =  ARGEMC 2, 2) 

C 

ARGEMC 2,1)  =  ARGEMC 1,2) 

ARGEMC 3,1)  =  ARGEMC 1,3) 

ARGEM(3,2)  =  ARGEM(2,3) 

ARGEMC4.1)  =  ARGEMC 1,4) 

ARGEMC4.2)  =  ARGEMC 2, 4) 
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ARGEM(4,3)  =  ARGEM(3,4) 

RETURN 

END 


SUBROUTINE  ASSEMB 


THIS  SUBROUTINE  ASSEMBLES  THE  ELEMENT  STIFFNESS,  ELEMENT  GEOM. 
STIFFNESS  MATRIX  TO  OBTAIN  THE  GLOBAL  MATRIX 


IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/ELEMNT/  EK(16,16,6),  EKT(16,16),  EG(16,I6) 

COMMON /ELEMNB/  ARSTF(4,4),  ARGEM(4,4) 

COMMON /PROPT1/  El(15),  E2(15),  G12(15),  P012(1S),  ANGL( 15) ,  Z(16) 

*  ,  NOD( 150,4),  NX,  NY,  NSDF,  LSDF(ISO) 

COMMON/PROPT3/  EBEAM,  ABEAM,  BMINER,  EPLT,  APLT,  THICK,  N0DB(18,2) 
COMMON/TOTAL  /  GK( 150 , 150 )  ,GG( 150 , 150 ) 

COMMON/BAND/  BANDK(45 , 150) ,  BANDKR(45 , 150) ,  BANDG(45 , 150) 
COMMON/EIG/NBAND ,  NREQ,  NEQ 
DIMENSION  NK(4) 

***  CLEAR  ARRAY 

DO  50  1=1,  NEQ 
DO  50  J=l,  NEQ 
GK( J , I )  =  0 . DO 
GG(J.I)  =  0 .DO 
50  CONTINUE 

NEL  =  NX*NY 
DO  100  N=1 ,  NEL 

DO  100  1=1,  4 

NR=(N0D(N,I)-1)*4 
DO  100  11=1,  4 
NR=NR+1 
L=( I- 1 )*4+II 
DO  200  J=l,  4 

NCL=(N0D(N,J)-1)*4 
DO  200  JJ=1,  4 
M=(J-1)*4+JJ 
NC=NCL+JJ 

GK(NR ,NC)=GK(NR ,NC)+EKT(L,M) 

GG(NR,NC)=GG(NR ,NC)+EG(L,M) 

200  CONTINUE 

100  CONTINUE 

***  IMPOSITION  OF  THE  BEAM  ELEMENT  MATRIX  *** 

DO  250  1=1,  NX 

NODI  =  4*(NODB( I , 1 ) - 1 ) 
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N0D2  =  4*(N0DB(I,2)-1) 

NK( 1)  =  NODI  +  1 

NK(2)  =  NODI  +  3 

NK(3)  =  NOD2  +  1 

NK(4)  =  NOD2  +  3 

DO  240  J=l,  4 
DO  230  K=1 ,  4 

GK(NK(K) ,NK( J) )  =  GK(NK(K) ,NK( J) )  +  ARSTF(K.J) 
GG(NK(K) ,NK(J) )  =  GG(NK(K) ,NK(J) )  +  ARGEM(K, J) 
230  CONTINUE 

240  CONTINUE 

250  CONTINUE 
C 

C  ***  IMPOSITION  OF  THE  BOUNDARY  CONDITIONS  *** 

C 

DO  310  NS=1,NSDF 
K=LSDF(NS)-NS+1 
DO  320  1=1, NEQ 
DO  320  J=K,NEQ 
GK( I , J)=GK( I , J+l) 

GG ( I , J ) =GG ( I , J + 1 ) 

320  CONTINUE 

C 

DO  330  I=K,NEQ 
DO  330  J=1 ,NEQ 
GK( I , J)=GK( 1+1 , J) 

GG( I , J)=GG( 1+1 , J) 

330  CONTINUE 

310  CONTINUE 
C 

C  ***  BAND  STORAGE  *** 

C 

NHBD  =  4*(NX  +  3) 

IF(NHBD.GE.NREQ)  NHBD  =  NREQ 
NBAND  =  NHBD  -  1 
C 

DO  20  J  =  1,  NREQ 

II  =  MAX0(1,  J-NBAND) 

DO  10  I  =  II,  J 
K  =  I-J+NBAND+1 
BANDK(K.J)  =  GK( I , J) 

BANDKR(K.J)  =  GK(I,J) 

BANDG(K , J)  =  GG( I , J) 

10  CONTINUE 
20  CONTINUE 
C 

RETURN 

END 

C 

C  _ 

C 

SUBROUTINE  EIGEN 

C  _ 

C 

C 

IMPLICIT  REAL*8(A-H ,0-Z) 
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COMMON /PR0PT1/  El(15),  E2(15),  G12( 15 ) ,  P012(15),  ANGL( 15) ,  Z(16) 

*  ,  N0D( 150,4) ,  NX,  NY,  NSDF ,  LSDF(150) 

COMMON/ PR0PT3/  EBEAM,  ABEAM,  BMINER,  EPLT,  APLT,  THICK,  N0DB(18,2) 
COMMON/ PR0PT4/  HT(15),  AAA,  BBB ,  NLYR,  NLYR1,  NLYR2 ,  NLYR3 
COMMON/TOTAL  /  GK( ISO , ISO) ,GG( 150 , 150) 

COMMON/BAND/  BANDK(45 , 150) ,  BANDKR(45 , 150) ,  BANDG(45 , 150) 

COMMON/E IG/NBAND,  NREQ,  NEQ 

COMMON/EIG2/  EGVC(150,5),  EGVC0( 150 ,5) ,  EGVL(5,4),  EGVL0(5,4), 

*  DTHK(50 ,5) ,  DGD(5) ,  EGVCB( 150 ,5) ,  IDMOD 
COMMON/NORMAL/  VNORMAL 

DIMENSION  WORK( 22500 ) ,  IND(3) 

EQUIVALENCE  (GG(1,1),  WORK(l)) 

EXTERNAL  OP,  IOVECT 

***  FACTOR  THE  BANDED  STIFFNESS  MATRIX  *** 

(LINPACK  SUBROUTINE  "DPBFA") 

LDA  =  45 

CALL  DPBFA(BANDK ,  LDA,  NREQ,  NBAND ,  INFO) 

IF( INFO. NE . 0)  STOP 

***  INPUT  FOR  SUB-SPACE  ITERATION  *** 

NVAL  =  3 

NFIG  =  10 
NPERM  =  0 
NMVAL  =  3 
NMVEC  =  150 
NBLOCK  =  NREQ/ 6 
IFCNBLOCK.GE.8)  NBLOCK  =  8 
MAXOP  =  NREQ 
MAXJ  =  150 

NREQ 2  =  2*NREQ 

***  CLEAR  THE  WORK  SPACE 

NV  =  IABS(NVAL) 

NWORK1  =  2*NREQ*NBLOCK  +  MAXJ*(NBLOCK+NV+2)  +  2*NBLOCK**2  +3*NV 
NWORK2  =  N*NBLOCK 

NWORK3  =  MAXJ*(2*NBLOCK+3)  +  2*NV  +  6  +  (2*NBLOCK+2)*(NBLOCK+l) 
IF(NWORK2 .GT .NWORK3)  NWORK3  =  NWORK2 
NWORK  =  NWORK1  +  NWORK3 

-  INITIAL  VALUE  OF  THE  EIGENVECTORS 

NPR7  =  NV 

IF(NBLOCK.LE.NV)  NPR7  =  NBLOCK 
DO  608  1=1,  NPR7 
DO  607  J=l,  NREQ 

WORKC (I-1)*NREQ+J)  =  EGVCB(J,I) 

607  CONTINUE 

608  CONTINUE 
C 

DO  1  1=1,  NREQ 

WORK(I)  =  0.D-05 
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CONTINUE 


1 
C 

CALL  DNLASO(OP,  IOVECT,  NREQ ,  NVAL,  NFIG ,  NPERM , 

*  NMVAL,  EGVL,  NMVEC ,  EGVC,  NBLOCK ,  MAXOP,  MAXJ , 

*  WORK,  IND,  IERR) 

C 

PRINT  1400,  IERR 

1400  FORMAT!  '  ERROR  PARAMETER  =',  15 ) 

PRINT  1420,  IND(l) 

1420  FORMAT( '  NO.  OF  CALLS  OF  OP  = ' , 15 ) 

PRINT  1430,  NPERM 

1430  FORMAT! '  NO.  OF  KNOWN  EIGEN  VALUES  = ’ ,15) 

IF(IERR.GT.O)  STOP 
C 

WRITE(6 ,662) 

662  FORMAT!//, 15X, 'EIGEN  SOLUTION  FOR  THE  PLATE', 

*  /  15X ,  35('=’)  ) 

C 

C  ***  POSTPROCESS  OF  THE  EIGENVALUES  *** 

C 

DENO  =  EPLT*APLT  +  EBEAM*ABEAM 
DO  760  1=1,  NV 

EGVL! 1,1)  =  1 . DO/EGVL! 1,1) 

SIG1  =  EGVL( I ,  l)*EPLT/DENO 
SIG2  =  EGVL( I , l)*E_EAM/DENO 
PI  =  SIG1*APLT 
P2  =  SIG2*ABEAM 
VNX  =  SIG1*THICK 
C 

DO  250  1=1,  NV 
C  ---  SAVE  EGVCB 

DO  245  K=1 ,  NREQ 

EGVCB!K,I)  =  EGVC(K ,  I ) 

245  CONTINUE 

C 

C  ***  OBTAIN  X  USING  THE  SUBROUTINE  "DPBSL2"  *** 

C 

CALL  DPBSL2 ! BANDK , 45 , NREQ , NB AND , EGVC ( 1 . I) ) 

C 

EGVL! 1,1)  =  EGVL! I , 1 )*VNORMAL 
WRITE! 6 , 666)  I,  EGVL(I.l) 

250  CONTINUE 
RETURN 
END 
C 

n 

w 

C  _ 

C 

SUBROUTINE  G?!N,M,P,Q) 

C  _ 

c 

IMPLICIT  REAL*8( A-H.O-Z) 

DIMENSION  P!N,M),  Q!N,M) 

DIMENSION  PBAR! 150 ) 

COMMON/BAND/  BANDK( 45 , 150 ) ,  BANDKR(45  150),  BANDG(45 , 150) 
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COMMON/EIG/NBAND,  NREQ ,  NEQ 
C 

DO  999  L=l,  M 

DO  105  J=l,  N 

PBAR(J)  =  P(J,L) 

105  CONTINUE 

C 

C  ***  OBTAIN  PBAR(N ,M)  USING  THE  SUBROUTINE  "DPBSL2"  *** 
C 

CALL  DPBSL2(BANDK ,45 ,NREQ ,NBAND , PBAR( 1) ) 

C 

C  ***  BANDED  GEOMETRIC  STIFFNESS  MATRIX  TIMES  PBAR  *** 

C 

NBAND1  =  NBAND  +  1 

LC  =  45 

CALL  BNDMUL(BANDG , PBAR , NBAND 1 , NREQ ,Q( 1 ,L) ,LC) 

C 

C  ***  OBTAIN  Q(N,M)  USING  THE  SUBROUTINE  "DPBSLl”  *** 

C 

CALL  DPBSLKBANDK, 45, NREQ, NBAND ,Q(1,L)) 

999  CONTINUE 
RETURN 
END 
C 

C  _ _ _ _ 

c 

SUBROUTINE  IOVECT(N ,M ,Q , J ,K) 

C _ 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/TOTAL  /  GK( 150 , 150) ,GC( 150 , 150) 

DIMENSION  Q(N,M) 

IF(K.EQ.l)  GO  TO  30 
DO  20  L=1 ,  M 
LI  =  J-M+L 
DO  10  1=1,  N 

GK(I,L1)  =  Q( I ,L) 

10  CONTINUE 

20  CONTINUE 

RETURN 
C 

30  DO  50  L=1 ,  M 

LI  =  J-M+L 
DO  40  1=1,  N 

Q( I ,L)  =  GK(I.Ll) 

40  CONTINUE 

50  CONTINUE 

RETURN 
END 
C 

C  _ _ _ — - 

c 

SUBROUTINE  BNDMUL( BAND ,WK .NBANDl , NREQ ,WK2 ,LC) 

C _ _ _ _ _ _ — - 

c 

IMPLICIT  REAL*8(A-H,0-Z) 
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DOUBLE  PRECISION  BAND(LC ,*) ,  WK(*) ,  WK2(*) 
NBAND  =  NBAND1  -  1 
DO  4200  1=1,  NREQ 
SUM  =  0.D0 
DO  4150  J=l,  NBAND1 
IPR  =  I  -  NBAND1  +  J 
IF(IPR.LE.O)  GO  TO  4150 
SUM  =  SUM  +  BAND( J , I )*WK( IPR) 

4150  CONTINUE 

DO  4160  K=1 ,  NBAND 
IPR  =  I  +  K 

IF( IPR. GT. NREQ)  GO  TO  4160 
SUM  =  SUM  +  BAND(NBAND1-K , IPR)*WK( IPR) 
4160  CONTINUE 

WK2( I )  =  SUM 
4200  CONTINUE 
RETURN 
END 


SUBROUTINE  GRMESH(NNM ,NEM) 


THIS  SOUBROUTINE  GENERATES  THE  BOOLEAN  CONNECTIVITY  MATRIX 


IMPLICIT  REAL*8(A-H,0-Z) 

COMMON /PROPT1/  El(15),  E2(15),  G12(15),  P012(15),  ANGLC15),  Z(16) 
*  ,  NOD( 150,4),  NX,  NY,  NSDF,  LSDF(150) 

C 

NX1  =  NX  +1 
NY1  =  NY  +1 
NXX1=NXX+1 
NYY1=NYY+1 
NEM  =  NX*NY 
NNM  =  NX1*NY1 
C 

NOD(l.l)  =  1 

NOD( 1,2)  =  2 

NOD( 1,3)  =  NX1  +  2 

NOD( 1 ,4)  =  NX1  +  1 

IF(NY.EQ. 1)  GO  TO  230 

M=1 

DO  220  N=2 ,  NY 

L  =  (N-1)*NX  +  1 
DO  210  1=1,  4 

NOD(L.I)  =  NOD(M , I )  +  NX1 
210  CONTINUE 

M=L 

220  CONTINUE 
230  IF(NX.EQ.l)  GO  TO  270 
DO  260  NI=2 ,  NX 
DO  240  1=1,  4 
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NOD(NI,I)  =  NOD(NI-l.I)  +  1 
240  CONTINUE 
M=NI 

DO  260  NJ=2 ,  NY 

L  =  (NJ-1)*NX  +  NI 
DO  250  J=l,  4 

NOD(L.J)  =  NOD(M ,  J)  +  NX1 
250  CONTINUE 

M=L 

260  CONTINUE 
270  RETURN 
END 
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APPENDIX  B.  OPTIMIZATION  PROBLEM  CODE  LISTING  USING 

IMSL 


PROGRAM  MAINIM 


VARIABLE  DECLARATIONS 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLEPRECISION  XLB(18),  XUB(18),  XGUESS(18),  T(18) 
PARAMETER  (IBTYPE=0,  IPRINT=3,  M=5 ,  MAXITN=100,  ME=1) 
EXTERNAL  FCN ,  GRAD,  INPUTB,  LIMITS,  GUESS,  DNCONG 


COMMON  BLOCKS 


COMMON/ELEMNT/ 
COMMON/ELEMNB/ 
COMMON/PROPTO/ 
COMMON /PR0PT1/ 

2 

COMMON/ PROPT2/ 

COMMON/ PR0PT3/ 

COMMON/ PR0PT4/ 

COMMON/TOTAL/ 

COMMON/BAND/ 

COMMON/EIG/ 

COMMON/EIG2/ 

2 

COMMON/FSAVE/ 

COMMON/BEAM/ 

COMMON/BEAM2/ 

COMMON/BIMO/ 

COMMON /LEVY/ 

COMMON/NORMAL/ 


EK( 16,16,6),  EKT( 16,16),  EG(16,16) 

ARSTF(4,4) ,  ARGEM(4,4) 

THETA,  ALEN,  BLEN 

E1C 15) ,  E2( 15)  ,  G12( 15)  ,  P012(15),  ANGL( 15 ) ,  Z(16), 
NOD( 150 ,4) ,  NX,  NY,  NSDF,  LSDF(150) 

DM( 3,3) 

EBEAM,  ABEAM,  BMINER ,  EPLT,  APLT,  THICK,  NODB(18,2) 
HT( 15) ,  AAA,  BBB,  NLYR,  NLYR1 ,  NLYR2 ,  NLYR3 
GK( 150 , 150) ,  GGC150.150) 

BANDK(45 , 150) ,  BANDKR(45 , 150) ,  BANDG(45 , 150 ) 

NBAND ,  NREQ,  NEQ 

EGVC( 150,5) ,  EGVCO( 150 , 5 ) ,  EGVL(5,4),  EGVLO(5,4), 
DTHK(50 ,5) ,  DGD(5) ,  EGVCB( 150 ,5)  ,  IDMOD 
VVV2( 100) 

E1B ,  E2B ,  G12B ,  P012B,  BBM,  HBM 

DBTHK(SO),  PBEAM 

VECRAT(5 ,5) ,  V99,  INDVC(5) 

SD11,  SD12,  SD22,  SD66,  ALPM,  ALPMS,  PEM,  OMEGA 
VNORMAL 


VNORMAL  IS  SCALING  FACTOR  FOR  T(NLYR3)  TO  INCREASE  ACCURACY 


VNORMAL  =  1.0D  03 


CALL  INPUTB 

N  =  NLYR  +  3 

CALL  LIMITS  (XLB ,  XUB) 

CALL  GUESS  (XGUESS) 

CALL  DNCONG  (FCN,  GRAD,  M,  ME,  N,  XGUESS,  IBTYPE,  XLB,  XUB, 
2  IPRINT,  MAXITN ,  T,  PMAX) 

END 
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SUBROUTINE  FCN  (M,  ME,  N,  T,  ACTIVE,  P,  G) 


IMPLICIT  REAL*8(A-H,0-Z) 

LOGICAL  ACTIVE(*) 
DOUBLEPRECISION  G(*),  T(*) 
EXTERNAL  INPUTC,  EIGANL4,  LEVYS 


COMMON  BLOCKS 


COMMON/ELEMNT/ 

COMMON/ELEMNB/ 

COMMON/PROPTO/ 

C0MM0N/PR0PT1/ 

2 

C0MM0N/PR0PT2/ 

C0MM0N/PR0PT3/ 

COMMON/ PR0P74/ 

COMMON/TOTAL/ 

COMMON/BAND/ 

COMMON/EIG/ 

C0MM0N/EIG2/ 

2 

COMMON/BEAM/ 

C0MM0N/BEAM2/ 

UPDATE  VALUES  OF  HT(*) ,  BBM ,  AND  HBM 

DO  100  1=1 ,NLYR 
HT( I )  =  T( I ) 

100  CONTINUE 

BBM  =  T(NLYRl) 

HBM  =  T(NLYR2) 

CALCULATE  VALUES  OF  PI,  P2,  P3 ,  AND  PBEAM  FOR  CURRENT  VALUES  OF 
HT(*) ,  BBM,  AND  HBM 

CALL  INPUTC 
CALL  EIGANL4 
CALL  LEVYS 
PI  =  EGVL( 1,1) 

P2  =  EGVL(2 , 1) 

P3  =  EGVL(3, 1) 

FUNCTION  DEFINITION 


EK( 16,16,6),  EKT( 16,16),  EG(16,16) 

ARSTF(4 ,4) ,  ARGEM(4,4) 

THETA,  ALEN,  BLEN 

El(  15 ) ,  E2( 15 )  ,  G12( 15 )  ,  P012(15),  ANGL( 15)  ,  Z(16), 
N0D( 150,4),  NX,  NY,  NSDF,  LSDF(150) 

DM( 3 ,3) 

E~EAM,  ABEAM,  BMINER ,  EPLT ,  APLT,  THICK,  N0DB(18,2) 
<  15),  AAA,  BBB,  NLYR ,  NLYR1 ,  NLYR2 ,  NLYR3 

GK( 150 , 150) ,  GG( 150 , 150) 

BANDK(45 , 150) ,  BANDKR(45 , 150)  ,  BANDG(45 , 150) 

NBAND,  NREQ,  NEQ 

EGVC( 150 ,5) ,  EGVCO( 150 ,5) ,  EGVL(5,4),  EGVL0(5,4), 
DTHK(50 , 5 ) ,  DGD(5) ,  EGVCB( 150 , 5 ) ,  IDMOD 
E1B ,  E2B ,  G12B,  P012B,  BBM,  HBM 
DBTHK(50) ,  PBEAM 


P  =  -T(NLYR3) 

CALCULATE  CONSTRAINT  VALUES  FOR  CURRENT  VALUES  OF  HT(*) ,  BBM,  &  HBM 

IF  (ACTIVE! 1))  G( 1)  =  APLT  +  ABEAM  -  THETA 
IF  (ACTIVE! 2) )  G( 2)  =  PI  -  T(NLYR3) 

IF  (ACTIVE(3) )  G(3)  =  ,999*P2  -  T(NLYR3) 

IF  (ACTIVE(4))  G(4)  =  ,998*P3  -  T(NLYR3) 
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IF  (ACTIVE(5))  G(5)  =  PBEAM  -  T(NLYR3) 
C 

RETURN 


SUBROUTINE  GUESS  (XGUESS) 


IMPLICIT  REAL*8(A-H ,0-Z) 

DOUBLEPRECISION  XGUESS(18) 

COMMON/ PR0PT4/  HT(15),  AAA,  BBB ,  NLYR ,  NLYR1 ,  NLYR2 ,  NLYR3 
COMMON/NORMAL/  VNORMAL 

READ  IN  XGUESS  VALUES 

DO  200  1=1 ,NLYR3 

READ(8,*)  XGUESS ( I ) 

0  CONTINUE 

SCALE  T(NLYR3)  GUESS  FOR  ADDITIONAL  ACCURACY 
XGUESS (NLYR3)=XGUESS(NLYR3)*VN0RMAL 
RETURN 


SUBROUTINE  LIMITS  (XLB ,  XUB) 


IMPLICIT  REAL*8( A-H ,0-Z) 

COMMON/PROPTO/  THETA,  ALEN,  BLEN 

COMMON/ PR0PT4/  HT(15),  AAA,  BBB,  NLYR,  NLYR1 ,  NLYR 2 ,  NLYR3 
DOUBLEPRECISION  XLB(18),  XUB(18) 

SET  LOWER  LIMITS  FOR  DESIGN  VARIABLES,  T(*) 

DO  100  1=1, NLYR 
XLB(I)  =  1.0D-03 
100  CONTINUE 

XLB ( NLYR 1)  =  1.0D-03 
XLB(NLYR2)  =  1.0D-03 
XLB(NLYR3)  =  1.0D-06 

SET  UPPER  LIMITS  FOR  DESIGN  VARIABLES,  HT(*) 

DO  200  1=1, NLYR 
XUB(I)  =  2.0D-01 
200  CONTINUE 

XUB<NLYR1)  =  2.0D-01 
XUB(NLYR2)  =  2.0D-01 
XUB(NLYR3)  =  1.0D  06 
C 

RETURN 
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APPENDIX  C.  OPTIMIZATION  PROBLEM  CODE  LISTING  USING 

DOT 


PROGRAM  MNDOT 


VARIABLE  DECLARATIONS 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLEPRECISION  XLB(18),  XUB(18),  XGUESS(18),  T(18),  DG(6,18), 
2  G(6) ,  AAC18.6),  WK(2000),  RPRM(20) 

DIMENSION  IPRM( 20) ,  IWK(IOOO) 

EXTERNAL  FCNDOT ,  INPUTB,  LIMITS,  GUESS,  DOT 


COMMON  BLOCKS 


COMMON /ELEMNT/ 
COMMON / ELEMN  B / 
COMMON /PROPTO/ 
COMMON / PROPT 1 / 

2 

COMMON/ PROPT2/ 
COMMON/ PROPT3/ 
COMMON/ PR0PT4/ 
COMMON/TOTAL/ 
COMMON/ BAND/ 
COMMON/E IG/ 
COMMON /EIG 2/ 

2 

COMMON/FSAVE/ 
COMMON/BEAM/ 
COMMON/ BEAM2/ 
COMMON/ BIMO/ 
COMMON/LEVY/ 
COMMON/NORMAL/ 


EK( 16,16,6),  EKT( 16 , 16) ,  EG(16,16) 

ARSTF(4 ,4) ,  ARGEM(4 ,4) 

THETA,  ALEN,  BLEN 

El( 15 ) ,  E2( 15 )  ,  G12( 15 ) ,  P012(15),  ANGL( 15 ) ,  Z(16), 
NOD( 150 ,4) ,  NX,  NY,  NSDF,  LSDFC 150) 

DM ( 3 , 3 ) 

EBEAM,  ABEAM,  BMINER,  EPLT,  APLT ,  THICK,  NODB(18,2) 
HT( 15) ,  AAA,  BBB,  NLYR,  NLYR1,  NLYR2 ,  NLYR3 
GK( 150 , 150) ,  GG( 150 , 150) 

BANDK(45 , 150 ) ,  BANDKR(45 , 150) ,  BANDG(45 , 150) 

NBAND,  NREQ,  NEQ 

EGVC( 150,5) ,  EGVCO( 150 ,5) ,  EGVL(5,4),  EGVLO(5,4), 
DTHK(50 ,5) ,  DGD(5) ,  EGVCBC 150 ,5)  ,  IDMOD 
VVV2( 100) 

E1B ,  E2B,  G12B,  P012B,  BBM,  HBM 

DBTHK(5Q) ,  PBEAM 

VECRAT(5 ,5) ,  V99,  INDVC(5) 

SD11,  SD12,  SD22,  SD66,  ALPM,  ALPMS,  PEM,  OMEGA 
VNORMAL,  SCALTH 


VNORMAL  IS  SCALING  FACTOR  FOR  T(NLYR3)  TO  INCREASE  ACCURACY 


VNORMAL  =  1.0D+02 


DEFINE  WORK  SPACE  REQUIREMENTS 

NRWK  =2000 
NRIWK  =  1000 

ZERO  RPRM  AND  IPRM 

DO  50  1=1,20 

RPRM(I)  =  0.0D0 
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IPRM(I)  =  0 
50  CONTINUE 

READ  IN  PROBLEM  DATA  AND  VALUES  OF  T(*),  XLB(*),  AND  XUB(*) 

CALL  INPUTB 

CALL  LIMITS  (XLB ,  XUB) 

CALL  GUESS  (XGUESS) 

DO  100  1=1 ,NLYR3 
T ( I )  =  XGUESS(I) 

100  CONTINUE 

SPECIFY  THAT  GRADIENTS  ARE  TO  BE  PROVIDED 
SET  PARAMETER  IPRM(l)  =  1 

IPRM(l)  =  1 

DEFINE  METHOD,  NDV,  NCON,  IPRINT,  MINMAX 
METHOD  =  1  FOR  MFD 
METHOD  =  2  FOR  SLP 

METHOD  = 

NDV  =  NLYR3 
NCON  =  5 
IPRINT  =  5 
MINMAX  =  0 

READY  TO  OPTIMIZE 

INFO  =  0 

200  CALL  DOT  ( INFO .METHOD , IPRINT, NDV , NCON ,T, XLB , XUB ,P, MINMAX ,G , 

2  RPRM , IPRM.WK ,NRWK , IWK ,NRIWK) 

PROVIDE  GRADIENTS  IF  DOT  REQUESTS  THEM 

IF( INFO .EQ .2)  GO  TO  300 

EXIT  IF  CONVERGENCE  IS  OBTAINED 

IF(INFO.EQ.O)  GO  TO  1000 

EVALUATE  OBJECTIVE  FUNCTION  AND  CONSTRAINTS 

CALL  FCNDOT  (T,  P,  G) 

GO  TO  200 

GRADIENT  EVALUATION 
300  CONTINUE 

CALCULATE  FUNCTION  GRADIENT  WRT  DESIGN  VARIABLES  AT  POINT  HT(*) 

DO  400  1=1 ,NI YR2 
WK( I)  =  0.0D0 
400  CONTINUE 
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WKCNLYR3)  =  -1.0D0 


IF  NUMBER  OF  ACTIVE  CONSTRAINTS  EQUAL  ZERO  RETURN  TO  DOT 

NGT  =  IPRM(20) 

IF(NGT.EQ.O)  GO  TO  200 

CALCULATE  CONSTRAINT  GRADIENTS  WRT  DESIGN  VARIABLES  AT  POINT  HT(*) 

DO  510  1  =  1  ,NLYR2 

DG( 1,1)  =  -DTHK( 1,1) 

510  CONTINUE 

DG( 1 ,NLYR3)  =  1 . 0D0 

DO  520  1=1 ,NLYR2 

DG(2,I)  =  -DTHK( I ,2)  *  0.999D0 
520  CONTINUE 

DG(2,NLYR3)  =  1.0D0 

DO  530  1=1 ,NLYR2 

DG(3 , I )  =  -DTHK(I,3)  *  0.998D0 
530  CONTINUE 

DG(3,NLYR3)  =  1.0D0 

DO  540  1=1 ,NLYR2 

DG(4 , I )  =  -DBTHK(I) 

540  CONTINUE 

DG(4 ,NLYR3)  =  1.0D0 

DO  550  1=1 ,NLYR 

DG(5 , I )  =  2.0D0  *  BLEN/THETA 
550  CONTINUE 

DG(5 ,NLYR1)  =  2.0D0  *  T(NLYR2) /THETA 
DG(5 ,NLYR2)  =  2.OD0  *  T(NLYRl) /THETA 
DG(5  ,NLYR3)  =  O.ODO 

STORE  ACTIVE  GRADIENT  INFORMATION  IN  ARRAY  WK(*) 

DO  600  K=1 ,NGT 
J  =  IWK(K) 

AA( 1 ,K)  =  DG( J , 1) 

AA(2 ,K)  =  DG(J ,2) 

AA(3,K)  =  DG(J,3) 

AA(4,K)  =  DG( J ,4) 

AA( 5 ,K)  =  DG( J,5) 

600  CONTINUE 

N1  =  NDV 

DO  800  K=1 ,NGT 
DO  700  1=1, NDV 

WK( I+Nl)  =  AA( I ,K) 

700  CONTINUE 

N1  =  N1  +  NDV 
800  CONTINUE 

GO  TO  200 
1000  CONTINUE 
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C 

C 


END 


SUBROUTINE  FCNDOT  (T,  P,  G) 


IMPLICIT  REAL*8(A-H,0-Z) 

DOUBLEPRECISION  G(*),  T(*) 

EXTERNAL  INPUTC,  EIGANL4,  LEVYS 

COMMON  BLOCKS  *** ** *****  *  ***** * ****** a-* 

COMMON/ELEMNT/  EK(16,16,6),  EKT(16,16),  EG(16,16) 

COMMON/ELEMNB/  ARSTF(4,4),  ARGEM(4 ,4) 

COMMON /PROPTO/  THETA,  ALEN,  BLEN 

COMMON / PROPT 1 /  El(15),  E2(15),  G12(15),  P012(15),  ANGL( 15) ,  Z(16), 
2  NOD( 150 ,4) ,  NX,  NY,  NSDF,  LSDF(ISO) 

COMMON / PROPT2 /  DM(3,3) 

COMMON / PR0PT3 /  EBEAM,  ABEAM,  BMINER ,  EPLT,  APLT,  THICK,  NODB(18,2) 
C0MM0N/PR0PT4/  HT(15),  AAA,  BBB,  NLYR ,  NLYR1,  NLYR2 ,  NLYR3 
COMMON/TOTAL/  GK(150,150),  GG(150,150) 

COMMON/BAND/  BANDK(45 , 150) ,  BANDKR(45 , 150) ,  BANDG(45 , 150) 
COMMON/EIG/  NBAND ,  NREQ ,  NEQ 

COMMON/EIG2/  EGVC(150,5),  EGVCO( 150 ,5) ,  EGVL(5,4),  EGVLO(5,4), 

2  DTHK(50 ,5) ,  DGD(5) ,  EGVCB( 150 ,5) ,  IDMOD 

COMMON/BEAM/  E1B ,  E2B,  G12B,  P012B,  BBM ,  HBM 
COMMON/BEAM2/  DBTHK(50) ,  PBEAM 
COMMON/NORMAL/  VNORMAL,  SCALTH 

UPDATE  VALUES  OF  HT(*) ,  BBM,  AND  HBM 

DO  100  1=1, NLYR 
HT( I )  =  T( I) 

100  CONTINUE 

BBM  =  T(NLYRl) 

HBM  =  T(NLYR2) 

CALCULATE  VALUES  OF  PI,  P2 ,  P3,  AND  PBEAM  FOR  CURRENT  VALUES  OF 
HT ( * ) ,  BBM,  AND  HBM 

CALL  INPUTC 
CALL  EIGANL4 
CALL  LEVYS 
PI  =  EGVL( 1,1) 

P2  =  EGVL(2 , 1) 

P3  =  EGVL(3 , 1) 

CALCULATE  FUNCTION  VALUE  FOR  CURRENT  VALUES  OF  HT(*) ,  BBM,  &  HBM 
P  =  -T(NLYR3) 

CALCULATE  CONSTRAINT  VALUES  FOR  CURRENT  VALUES  OF  HT(*) ,  BBM,  &  HBM 

G ( 1 )  =  (T(NLYR3)-P1) 

G(2)  =  (T(NLYR3) *0 . 999D0*P2) 

G(3)  =  (T(NLYR3)-0.998D0*P3) 
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G(4)  =  (T(NLYR3)-PBEAM) 

G(5)  =  (APLT  +  ABEAM) /THETA  -  1.0D0 
C 

RETURN 


SUBROUTINE  GUESS  (XGUESS) 


IMPLICIT  REAL*8( A-H ,0-Z) 

DOUBLEPRECISION  XGUESS(18) 

COMMON / PR0PT4/  HT(15),  AAA,  BBB,  NLYR ,  NLYR1 ,  NLYR2 ,  NLYR3 
COMMON/NORMAL/  VNORMAL 

READ  IN  XGUESS  VALUES 

DO  200  1=1 ,NLYR3 

READ(8 ,*)  XGUESS(I) 

200  CONTINUE 

SCALE  TCNLYR3)  GUESS  FOR  ADDITIONAL  ACCURACY 
XGUESS ( NLYR3 ) =XGUESS ( NLYR  3 ) *VNORMAL 
RETURN 


SUBROUTINE  LIMITS  (XLB ,  XUB) 


IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/ PROPTO/  THETA,  ALEN,  BLEN 

COMMON /PR0PT4/  HT(15),  AAA,  BBB,  NLYR,  NLYR 1 ,  NLYR 2 ,  NLYR3 
DOUBLEPRECISION  XLB(18),  XUB(18) 

SET  LOWER  LIMITS  FOR  DESIGN  VARIABLES,  T(*) 

DO  100  1=1, NLYR 
XLB(I)  =  1.0D-03 
100  CONTINUE 

XLB ( NLYR 1)  =  1.0D-03 
XLB(NLYR2)  =  1.0D-03 
XLB(NLYR3)  =  1.0D-06 

SET  UPPER  LIMITS  FOR  DESIGN  VARIABLES,  HT(*) 

DO  200  1=1, NLYR 
XUB(I)  =  2.0D-01 
200  CONTINUE 

XUB ( NLYR 1)  =  2.0D-01 
XUB(NLYR2)  =  2.0D-01 
XUB(NLYR3)  =  1.0D  06 
C 

RETURN 
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